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Dispersionless Single-Mode Lightguides 
With a Index Profiles 


By U. C. PAEK, G. E. PETERSON, and A. CARNEVALE 


(Manuscript received November 14, 1980) 


By means of a numerical solution to Maxwell’s equations, we 
calculate those parameters necessary to design and fabricate single- 
mode lightguides. These include optimum core radius and profile 
parameter a. In the design the dispersion is minimized by varying 
the core radius while the relative index difference A, the wavelength 
A, and profile parameter are held fixed. A comparison of calculated 
dispersion with experimental data shows excellent agreement. 


1. INTRODUCTION 


Recently, single-mode lightguides have been developed that achieve 
transmission losses as low as 0.5 dB/km and 0.2 dB/km at wavelengths 
of 1.3 wm and 1.55 pm, respectively. When the total lightguide disper- 
sion is reduced to zero at the operating wavelength, a transmission 
system can be realized with wide repeater spacings and extremely large 
bandwidths.”“* In fact, bandwidths in excess of 1 GHz/100 km are 
expected. Clearly, such a lightguide system would be ideal for undersea 
cable and other long-distance transmission applications. 

The design of single-mode lightguides with zero total dispersion 
requires an accurate description of this parameter in terms of profile 
shape, materials properties, and core radius. This necessitates a solu- 
tion of Maxwell’s equations for the fundamental HE), lightguide mode 
as well as the TE, and TMo; modes. 

For radially inhomogeneous media, it is usually not possible to 
obtain these solutions as analytical expressions of a closed form. To a 
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very large extent the design of single-mode lightguides has in the past’ 
focused on rectangular-shaped index profiles. This is because Max- 
well’s equations are very much easier to solve in this case. When 
nonrectangular profiles have been considered, numerous approxima- 
tions have often been employed. 

To avoid these difficulties, we have developed a numerical technique 
to obtain exact solutions to the vector form of Maxwell’s equations.” 
These equations are written as four coupled simultaneous first-order 
differential equations. The effective index N, for the single propagating 
HE, mode is found by solving the characteristic equation derived 
from the matching conditions of the field components at the core 
cladding interface. The dispersion is then calculated from N,. In 
addition, we find the conditions under which the TE; and TMo: modes 
propagate. 

Our computing procedures do not impose any restrictions on the 
index profile of the fiber.°® However, for the sake of simplicity and 
familiarity we choose the a-index profile. To confirm the analysis, we 
compare the computed results with known experiments when possible. 
The calculated results include the dispersion, the optimum core radius 
(defined as the value at which the zero total dispersion will occur for 
a given relative index difference, wavelength, and a value), and the 
corresponding cutoff frequency. 

To keep the accuracy as high as possible, material dispersion is 
included in the calculation from the onset. The subtle interaction 
between material properties and waveguide properties determine the 
propagation characteristics of the lightguide. Details of the numerical 
procedure can be found in the appendix. 


Il. DESIGN OF A DISPERSIONLESS SINGLE-MODE FIBER 


In the design of a single-mode fiber, the relative index difference A, 
the operating wavelength A, and the cladding index N2 are normally 
specified. The objective of the design is to make the total dispersion 
for the HE,, mode equal to zero, the total dispersion D; being defined 
as the time spread of a narrow pulse per fiber length L per source 
spectral width di. Qualitatively speaking, this can be accomplished by 
adjusting the radius and profile shape so that the waveguide dispersion 
exactly counterbalances the material dispersion. In addition, the design 
must be such that the TE, TM, and other higher-order HE and EH 
modes are not propagating. In the analysis shown in the appendix, the 
number of modes allowed to propagate in the fibert core are not 
restricted. However, the single-mode analysis is nothing but a partic- 
ular case of multimode propagation and it can be easily achieved 
simply by setting the angular mode number M = 0 and 1 in the 
elements of matrix A [eq. (18) in the appendix] while fixing the radial 
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mode number @ = 1. When M = 0, the cutoff frequency for the TE 
and TM can be determined. Below this cutoff point, only a funda- 
mental mode (HEji) propagate in the fiber. All the modes other than 
M = 0 are hybrid modes. Further details are explained below. 

When egs. (22) and (23) are substituted into eq. (21), determinant 
(21) decomposes into two uncoupled equations, giving the TM and TE 
modes, 

TM mode: 


(BiAa — K(fi)yoAn) _ 


alas Sita ET hg 1 
(Bi Az» = K(f) yoAr2) (1) 


TE mode: 
(BiAz2 x yoA42) | 
(BiAs: — yoAa) ; 


When Jj is the maximum value of the index of refraction in the core 
and Nz is the cladding index, the effective index of the bounded modes 
must satisfy the condition N2< N,. < N;. Therefore, from eqs. (1) and 
(2) the effective index N,(0, 1) for the TMo; and TEo: mode can be 
found.’ The notation N.(M, Q) refers to the effective index of a mode 
having angular mode number M and radial mode number Q. 

In practice, the core radius is changed until N.(0, 1) is equal to Ne. 
This gives us the cutoff radius a,- for the TEo; and TMo, modes. The 
normalized cutoff frequency is then calculated from 


Vox aris JN? — Ni. (3) 


(2) 





Now if V < V., only a single mode (HE};) propagates. The effective 
index N, for this mode is calculated in the single-mode region as a 
function of the radius and is the principal quantity used to determine 
the dispersion of the HE, mode. 

The relationship between the effective index N-(1, 1) and the group 
index N,j is given by 
dN. 
dh - 
The propagation time of a pulse through a fiber of length L and group 
index Nz is 





Nz =Ne— A (4) 


L 
t= G Ne (5) 


If there is a variation of N, with wavelength, there is a dispersion or 
spread in pulse width. Thus 

_LdNz 
~ C dx 





dt dn, (6) 
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where dA is the spectral width of the source. 
We can substitute N, from e (4) into eq. (6) and get 








d’N. 
dt = —-— A (dN) De’ (7) 
The total dispersion D; is 
dt 1 
D: = a . L (8) 
Thus, 
d\ a’Nn. 
Oa ane @) 


Obviously, we need to solve Maxwell’s equations for N. and then 
calculate d?N./dd? to obtain D,. 

Following the procedure described in our earlier papers, it is con- 
venient to describe the dispersive properties of the cladding by a 
modified Sellmeier formula.*” Thus, 


C3 2 C, ‘ Cs 
(2-1) (P-L)? Q?-D)” 


where / = 0.035. The coefficients C; are given in Table I. Figure 1 
shows the variation of the silica cladding’s refractive index, for the 
lightguides considered here. 

Since the relative index difference A is rather small, typically ranging 
from 0.2 percent to 0.8 percent, the core center index N; can be written 
as 


No = Co + Cid? + Cod4 + (10) 


= Ne 
N, = {a - Ay (11) 


To obtain a feel for the size of the numbers involved, the following 
data are useful: 


A = 1.33 um, 
Nez = 1.446925, 
N, = 1.449825 [as calculated by eq. (11)]. 


Table | 
Co =  1.4508554 
C, = —0.0031268 
C. = —0.0000381 
C; = 0.0030270 
C, = —0.0000779 
C; = 0.0000018 
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1.46 


1.45 


INDEX OF REFRACTION 


1.44 
0.6 0.8 1.0 1.2 1.4 1.6 


WAVELENGTH IN MICROMETERS 


Fig. 1—Refractive index of silica as a function of wavelength. 


Since we assume a wavelength dependence for Ne, described by eq. 
(10), Ni as calculated by eq. (11) is also wavelength dependent. The 
following profile formula is particularly useful and will be employed in 
the fiber design: 


N(r) = N,[1 - Ar*]. (12) 


A few comments seem in order: When a — ~ the profile is rectangular, 
when a = 2 it is parabolic, when a = 1 it becomes linear, and when 
a < 1 the profile develops a cusp. Figure 2 illustrates these shapes. 

In practice, one specifies A, A, and the profile parameter a. The 
computer then searches for a core radius that makes D; = 0 by means 
of Mueller’s iterative method. If it can be assumed that the total 
dispersion is the sum of waveguide dispersion D,,(A, a, a, A) and 
material dispersion D,,(A), then the physics of the problem is easy to 
understand. As will be shown later, this separation is quite accurate. 
Thus, a dispersionless single-mode fiber is one in which material 
dispersion is exactly balanced by waveguide dispersion. This is illus- 
trated by Fig. 3. 

Curve D,,,(A) shows the material dispersion as a function of wave- 
length for a single-mode lightguide with A = 0.2 percent and a rectan- 
gular profile. Note that it changes sign and passes through zero at 
about 1.27 ym. Also plotted are waveguide dispersion curves for three 
different radius values. Note that all three curves are negative and that 
as the radius decreases the magnitude of the dispersion increases at a 
given wavelength. . 

From Fig. 3 we can draw a number of conclusions. First, the 
lightguide with the materials properties shown can be made totally 
dispersion free only at wavelengths longer than 1.27 um. Second, as 
the wavelength gets longer the guide radius must get smaller. Finally, 
at longer wavelengths a much larger amount of material dispersion 
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0.002 


0.001 


N—No, REFRACTIVE INDEX DIFFERENCE 





0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 
R/la, NORMALIZED RADIAL COORDINATE 


cg ee a oe OE 


Fig. 2—Index profile shapes for various power laws. When a — © the profile is 
rectangular, when a = 2 it is parabolic, and when a = 1 it is linear. If a is less than 1, the 
profile develops a cusp. 


must be compensated for by waveguide dispersion. This requires 
greater precision in the waveguide parameters than when the guide is 
designed to operate at the zero of material dispersion. Hence, it is 
desirable to have a variety of materials available for lightguide design 
purposes. 

Finally, a typical total dispersion curve D,, illustrating lightguide 
operation at 1.3 zm, is plotted in Fig. 3. 


lll. SINGLE-MODE LIGHTGUIDES WITH a PROFILES FOR 1.33-um 
OPERATION 


Single-mode lightguides that are designed to operate at 1.33 um with 
A = 0.2 percent are of much current interest. The cutoff frequencies 
V.(a) for the a-index profiles can be found from eqs. (1) and (2). Figure 
4 plots the results. There is little change in V, until the profile 
parameter is below 10, then it increases rapidly. At a = 1, Ve and Vopr: 
are 4.383 and 2.724, respectively. The data of Fig. 4 agree well with 
that in the literature.®” 

To determine the radius a, for zero total dispersion, a computer 
search is done for values less than a., where a, is the cutoff radius [see 
eq. (3)]. Figure 5 shows the result for three a values, namely 100, 2, 
and 1. The radii turn out to be 4.142, 5.725, and 6.294, respectively. 
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Note that the total dispersion D; is most sensitive to the radius for 
a = 100, and least sensitive when a = 1. For design purposes, we plot 
@opt as a function of the profile parameter a in Fig. 6. This plot also 
includes the cutoff wavelength A., where A, is defined by the following 
formula: 


Ve = (Vopt/ Ve) “Xr. (13) 


Note that both V.,: and V, are employed in the calculation. Therefore, 
to allow only one mode to propagate in the core, the operating 
wavelength must be longer than the cutoff wavelength .. For a step 
index profile with A = 0.2 percent and A = 1.33 pm, A, is close to 1 ym. 


7 Dy \X,a,a@,A) 
a 


DISPERSION (ps/km - nm) 





1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 
WAVELENGTH IN MICROMETERS 


Fig. 3—A plot showing material dispersion D,,, waveguide dispersion D,,, and total 
dispersion D,. Note that because of the difference in sign between the material dispersion 
and the waveguide dispersion it is possible for one to cancel the other at a particular 
wavelength. 
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Vo (a), CUTOFF FREQUENCY 
Vopt( @), OPTIMUM FREQUENCY 





1 2 4 6 8 10 2 4 6 8100 2 4 6 
a, PROFILE EXPONENT 


Fig. 4—A plot of V. and Vo as a function of the profile parameter. 


The calculations described thus far in this section have made no 
assumption concerning the separation of the total dispersion D; into 
its component parts D,, and D,. Our numerical procedure simply 
searches for a zero in total dispersion and does not assume that D; = 
Dn + Dw. 

However, the material dispersion contribution can be calculated 
directly by means of eqs. (10) and (9), where Nz is substituted for N-. 
The waveguide dispersion can be found by solving Maxwell’s equations 
numerically via the procedure outlined in the appendix, with N2 being 
wavelength independent. 

Figures 7a through 7c show the results for a = 100, 2, and 1. In all 
three figures the material dispersion, which is of course independent 
of core radius, is represented by the horizontal dashed line. The 
negatives of the waveguide dispersions are dependent on core radius 
and are the solid lines in the figures. The intersection of the solid and 
dashed lines yields the optimum radii. As is evident from the figures, 
the radii are identical with the previous calculations. 


IV. A COMPARISON WITH EXPERIMENTAL DATA 


Miya, Terunuma, and Hosaka” have fabricated GeO» doped single- 
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mode lightguides for the 1.3 um and 1.5 pm region. Their experimental 
data are very useful to check the theory presented in this paper. They 
give data on a fiber with A = 0.2 percent and a zero total dispersion at 
A = 1.33 um. A second fiber has A = 0.74 percent and a zero total 
dispersion at \ = 1.54 pm. Both are step-index fibers. Their experi- 
mental data are measured with a fiber Raman laser excited by a Q 
switched, mode-locked Nd:YAG laser. The spectral range studied is 
1.1 to 1.7 pm. The optimum radius is 4.142 um for A = 0.2 percent and 


Dt, TOTAL DISPERSION IN ps/km-nm 





“0 2 4 6 8 10 12 
a, CORE RADIUS IN MICROMETERS 
Fig. 5—Total dispersion versus core radius for lightguides with a = 1, a = 2, and 


a= 100. Note that the total dispersion curve is steepest for a = 100 and shallowest for 
a=1. 
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0.8 


aopt. OPTIMUM CORE RADIUS 
IN MICROMETERS 
oO 
© 
Ac, CUTOFF WAVELENGTH IN MICROMETERS 





° 
N 


2 ee 
1 2 4 6 8 10 2 4 6 8100 2 4 6 
qa, PROFILE EXPONENT 


Fig. 6—A plot of cutoff wavelength A, and optimum radius aop: versus a. This data is 
very useful for design purposes. 


A = 1.33 ym. The optimum radius is 2.30 um for A = 0.75 percent and 
A = 1.53 ym. Figure 8 compares their experimental results (circles and 
triangles) with our calculations (solid and dashed lines). The agreement 
is excellent. 


V. DISCUSSION 


At a wavelength of 1.33 pm, dispersion-free lightguides can be 
designed with improved properties over the rectangular profile guides. 
Our calculated results of a., and cutoff frequency V.(a) both show a 
rapid increase as the profile parameter a drops below 10. For a 
triangular index profile (a = 1), the core size is over 50 percent larger 
than the size of a step-index core. Thus, due to ease of handling, 
coupling to the source, and splicing, a larger core size is desirable for 
the same value of A and A. 

As the profile parameter a drops below 10, Vo, also begins to 
increase substantially. It is significant that in the case of a = 1, Vopt ~ 
2.7, which is a substantial increase over the step-index case, Vopt ~ 
1,.8,1)12 

In addition, Fig. 5 shows that a lightguide with a triangular profile 
is less susceptible to variations in core radius than a step-index guide. | 
Therefore, because of the larger core size and less sensitivity in its 
tolerance, a triangular profile may be less difficult to fabricate than a 
rectangular profile. 
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APPENDIX 

To study the propagating modes in a fiber, we first assume that the 
cross section of the fiber is concentric in the core and cladding and 
that the core diameter 2a is much smaller than the outer diameter D 
(Fig. 9). Further, we assume that the core is a lossless, radially 
inhomogeneous medium with a scalar permittivity «. A relative per- 
mittivity (dielectric constant) is defined by k = €/é€o, where €o is the 
value in free space. However, as is well known in a nonconducting 
medium, the permittivity becomes equal to the square of the index of 
refraction N. For source-free fields in a dielectric waveguide, Maxwell’s 
equations reduce to 


oH 
VX E =p, 
bot 7” 
14 
a eee 
ot 


—Dy\A.a,a,A) 


DISPERSION IN ps/km.nm 





a, CORE RADIUS IN MICROMETERS 


Fig. 7(a)—A ay of waveguide dispersion and material dispersion as a function of 
core radius, a = 1 
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—Dy\A.a,a,A) 


DISPERSION IN ps/km-nm 





5 6 7 8 
a, CORE RADIUS IN MICROMETERS 


. Fig. 7(b)—A plot of waveguide dispersion and material dispersion as a function of 
core radius, a = 2. 


where E(R, ¢) and H(R, t) represent the electric and magnetic fields. 
lo is a scalar permeability in free space. 

In a cylindrical coordinate system R = {R, ¢, 2}, the solution of eq. 
(14) is described by the vector components of two fields, {Er, Ey, Ez} _ 
and {Hr, H;, Hz}. Among these components we are primarily inter- 
ested in finding the tangential components {F,, E.} and {H,, H-}. 
These will be continuous through the core-cladding interface. Once 
the tangential components are known, the radial components Er and 
Hr can be obtained from these components. To find a complete set of 
bounded modes we consider the following form of solution:”® 


E(p, $, 2, t) = E(p) exp[t(wt — Md — Bz)], (15) 


where w and M are angular frequency and mode number, and £ is 
propagation constant along the z axis. 
Introducing a new variable A, the field components are written as 
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Ai E, 


Ae —I1ZopH, 
A=} cea (16) 
Ag —1Z)Hz 


where p = KR = 27R/), Zo is the wave impedance defined by (j1o/€0),"” 
and J is a wavelength. 

Substituting eq. (16) into eq. (15) yields a set of the first-order 
ordinary differential equations.’** In the following equation, 


—— =-A(p)A, (17) 


the operator A(p) is a 4 X 4 matrix that contains important property 

constants and parameters, for instance, the index of refraction N, 

angular mode number M, and the effective index N, defined by B/k. 
The 4 x 4 matrix A(p) can be written as . 


—Dy\r,a,a,A) 


DISPERSION IN ps/km.nm 





a, CORE RADIUS IN MICROMETERS 


Fig. 7(c)—A plot of waveguide dispersion and material dispersion as a function of 
core radius, a = 1. 
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Dz, TOTAL DISPERSION IN ps/km-nm 





1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 
WAVELENGTH IN MICROMETERS 


Fig. 8—A comparison of theory (dashed and solid lines) with experiment (see text). 


0 (N2/x) —1 0 —MN./k 
p'k — M’ 0 MN. 0 18 
0 MN./« 0 0? — (M2/x)|° 8) 
—MN. 0 N2—« 0 


The method of numerical computation has been given in our earlier 
work,” which briefly gives the procedure for finding the solutions and | 
emphasizes what is needed to describe a single-mode fiber design. 

Equation (17) is a familiar type of differential equation. A fourth- 
order Runge-Kutta method is used to solve it numerically. A concise 
description is as follows. First, we know that there are two possible 
solutions in both the core and cladding, and they are linearly inde- 
pendent. Therefore, a general solution will be a linear combination of 
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two solutions, denoted by A; in the core and I; in the cladding. 
In the core region, 


A = Ay Ain + AoAie, (19) 
and in the cladding region, 
T; = A3Aj3 + AgAja. (20) 


The solutions A;; and Aj2 are transferred to the core-cladding inter- 
face through the operation of eq. (17), starting with two initial condi- 
tions given at the center of the core.*”* 

At the interface p = p;, each has to be matched with the solution in 
the cladding to satisfy the continuity condition mentioned earlier. 
Thus, the matching condition yields a set of homogeneous equations 
for the constants A;. Hence, the characteristic equation is the vanishing 
determinant of the system of equations, 


det| Aa | = 0. (21) 





Fig. 9—Cross section of a fiber and cylindrical coordinate eet R or p is the radial 
esoetinais. ¢ is the angle, and z is the axial direction (p = KR). 
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Particularly, considering M = 0, we can write the solutions in the core 
and cladding as follows. The two core solution vectors at p = p; are 
An and Aj,2(p;). Expressed in terms of their components, they are 


An Aye 
er ae Agi ee: Age 
An = An and Ai = ie (22) 
Ag Aae 
And the two cladding solutions at p = p; are taken as 
Bi 0 
| (Sa yo(Ss) | a 0 | 
Ais — 0 Wo(§:) and Au yo( i) Wo(6i), (23) 
0 Bi 


where 
Bi=[Ne—«(S)]'", — & = Bini, 
yo = Sie [Ko($i)/Ko(fi)], and Wo( Si) = Ko( Si) /Bi. 


The prime in the above equation denotes differentiation with respect 
to ¢. 
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Closed Markovian networks of queues that have the product form 
in their stationary probability distributions are useful in the perform- 
ance evaluation and design of computer and telecommunication 
systems. Therefore, the efficient computation of the partition func- 
tion—the key element of the solution in product form—has attracted 
considerable effort. We present a new and broadly applicable method 
for calculating the partition function. This method can be applied to 
very large networks, which were previously computationally intrac- 
table. Most of the paper details applications of this approach to a 
network class which arose in modeling an interactive processor. We 
show that the partition function and derivatives such as mean values 
(response times, cPU utilizations, etc.) may be represented by integrals 
and their ratios. The integrands contain a parameter N which is 
large for large networks. Next, the classical techniques of asymptotic 
analysis are applied to derive three main power series expansions in 
descending powers of N to correspond to normal, high, and very high 
usage. This work emphasizes multiple terms in the expansions for 
precision and error analyses. 


1. INTRODUCTION 


The theoretical results on the product form of the stationary distri- 
butions of large classes of Markovian queuing networks continue to 
have a profound influence on computer communications, computer 
systems analysis, and traffic theory.’* These results make at least 
feasible the analysis and synthesis of the large systems of ever increas- 
ing complexity being considered in these areas. The subclass of closed 
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networks of queues is more difficult to analyze than the open networks 
because there is no stationary independence of the network nodes. 
However, the incentive for investigating the closed networks does exist 
since they have been used to model multiple-resource computer sys- 
tems,” multiprogrammed computer systems,** time-sharing,” and 
window flow control in computer communication networks;*”° net- 
works with external inputs subject to blocking require the analysis of 
a large number of closed networks.’””” The closed network model that 
we shall use for illustrative purposes arose in the modeling of a central 
processor in a node of a computer network. This network is subject to 
a variety of processing demands. In recognition of the utility of closed 
networks, considerable research and commercial interest has been 
directed towards developing efficient procedures for computing the 
partition function (the normalizing constant), the only element of the 
product form solution requiring significant computation.'*"8 

However, as these existing recursive techniques are applied to the 
problems of particular interest in the Bell System, wherein the con- 
stituents of the closed chains are many and the number of chains are 
many, their shortcomings are observed to be severe in the amount of 
computing time and memory required and the accuracy attained. A 
more detailed account appears in Section 2.4 and Section IX. Briefly, 
the existing recursive techniques are largely ineffectual. 

We present a new way to view the problem, which surmounts many 
of the difficulties associated with large networks. The approach is 
broadly applicable—as indicated in Section X—even though the paper 
is a detailed account of applications to a specific class of closed 
networks. The new approach consists, first, of recognizing that the 
partition function may be written as an integral with a large parameter 
N present in the integrand to reflect the large size of the network. 
Next, the classical techniques of asymptotic analysis are applied to 
derive an asymptotic power series, typically in descending powers of 
N. The integrand will have fundamentally different properties in 
different ranges of the system parameters and this will require corre- 
spondingly different expansions. Thus, in this paper we develop three 
separate series expansions—Proposition 3 (Section IV), Proposition 12 
(Section VII), Proposition 17 (Section VIIIT)—each corresponding to a 
specific range of values of the usage parameter a. It is worth empha- 
sizing that, commensurate with an objective of providing solutions 
with any desired accuracy, we give procedures for generating multiple 
terms in the asymptotic expansions, not just the dominant term. In 
Section VIII, we unify the preceding results by giving a common 
expansion that holds uniformly in the system parameters. The uniform 
expansions introduce in a natural way the parabolic cylinder (or 
Weber) functions, a classical family of special functions with many 
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antecedents and ties with other special functions.’* Besides duplicating 
the specialized expansions derived earlier, the uniform expansion 
makes available for use the many well-documented and tested expan- 
sions that are known in connection with parabolic cylinder functions. 

Section IX describes a user-oriented software package that has been 
written in C-language to implement the approach developed here. We 
supply results obtained by the package on four test problems that 
arose in analyzing performance of a Bell System project. Also reported 
are the results of a comparison with a well known, commercially 
marketed package that obtains solutions recursively. Our. package is 
able to solve the large problems, which are well beyond the range of 
the other package, and, surprisingly, solve the small problems as well 
with errors that have small bounds. 

Section X provides the basis for extending the approach developed 
here to quite general multiprocessor, multidiscipline queuing networks. 
We show that for most networks that have been shown to have the 
product form in their stationary probability distribution, the partition 
function has an integral representation. The expansions appropriate 
for its computation are not considered here. 

Not surprisingly, the new representation of the partition function as 
an integral—the starting point of our computational procedures—may 
be exploited anew to derive analytical estimates and bounds of the 
quantities of interest, such as throughput, mean response time, etc. 
We demonstrate particularly in Section 5.3 that these formulas explic- 
itly exhibit the system parameters and as such are rather useful as 
design and synthesis aids. (The bounds are also useful as checks on 
the computational procedure.) Purely computational procedures by 
themselves do not yield this particular form of insight into system 
behavior. 

The asymptotic sequences used typically are power series in N™', 
where recall N is the generic large parameter.’* Thus, the number of 
terms required to achieve the desired accuracy decreases with increas- 
ing N. In contrast, with recursive solutions the computational com- 
plexity grows with the network size. Also, the asymptotic methods 
handle increased numbers of classes of constituents with little incre- 
mental difficulty, while the computational complexity in recursive 
methods grows geometrically. Thus, the contrasting techniques are 
not replacements for each other but complementary: loosely speaking, 
the recursions are most effective for smaller networks, while the 
asymptotic expansions are most effective for large networks. 

The contrasting behavior with respect to a large number of classes 
is of particular importance in computer communications where, as 
Reiser,’ Schwartz,’ and others have pointed out, traffic corresponding 
to each source-destination pair is treated as a separate class and the 
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network closure follows from the windows employed in flow control. 
Reiser has developed heuristics to cope with this situation. 

References 11, 15, and 16 also contain results pertaining to large 
networks. 


il. NETWORK MODEL AND KNOWN RESULTS 
2.1 Model 


In the model (see Fig. 1) each constituent, which may be thought of 
as a terminal or station, of the closed network spends alternating 
periods of time in the two nodes that constitute the network—the 
‘think’ node (also node 1) and the ‘cpu’ node (node 2). The think time 
in each cycle for each constituent is an independent random variable 
with an exponential distribution. The time spent in the cpu node 
depends on many factors since it is here that there occurs interaction 
between constituents being serviced. We stipulate that the cpu disci- 
pline is ‘processor sharing’* and that the desired service time (i.e., the 
time required to service the job if the entire cpU was dedicated to the 
job) is an independent exponentially distributed random variable.'* 
The ‘think’ node is thus an o-server center and in the terminology of 
Ref. 4, nodes 1 and 2 are respectively Type 3 and 2 centers. 

We stipulate that there is sufficient statistical inhomogeneity 
amongst the constituents to justify the existence of several, say p, 
classes of constituents, with class 1 having K; constituents, 1 <7 p. 
Statistical homogeneity applies within a class in that pj and p;2 will 
respectively denote the mean think time and the mean desired service 
time that are common to all in class 7. The variations among these 
mean values may be quite substantial. 

Our involvement with this model arose while modeling behavior of 
traffic through a processing node of a computer network. The number 
of classes of constituents is at least five, namely, time sharing; inquiry/ 
response and data-base query; batch and remote job entry; messages 
and broadcast; data entry/collect and screen type jobs. The mean 
values {p;;} are obtained from benchmark measurements. In another 
variation of this problem, a finer classification of constituents was 
considered. Our interest is in cases where the individual class popula- 
tions {K;} extend to several hundred, while the number of classes is of 
the order of ten. 


* In the processor-sharing discipline there is no overt queuing because all, say n, jobs 
present in the node simultaneously receive service at 1/n times the rate given to a single 
job by the processor. Thus, the rate given to a single specific constituent fluctuates with 
time. This discipline is the limiting case of the round robin discipline as the time 
quantum given to each job becomes arbitrarily small. 
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TERMINALS 











CLASS 7 
CENTRAL 
PROCESSOR 
SYSTEM 
CLASS p 


CPU CPU 
THINK THINK 
TIME 


|. RESPONSE -»| 


TIME 


NODE 1 NODE 2 
co - SERVER PROCESSOR - SHARING 


ne 


(b) 


Fig. 1—(a) There are p classes of constituents—shown as terminals—with K; constit- 
uents in class j. (b) Constituents spend alternate periods of time in the think node and 
processor-sharing CPU node. 


2.2 Product form solution 


If N,; is the stationary random variable denoting the number of 
constituents of class i in node J, and if 7 is the following stationary 
probability n(n, eee, Mp1; Miz, *°*, Np2) = Pr{ Niu = N11, ***, Noi = 
Npi; Nie = Mie, +++ , Np2 = Np2], then it is known that with the left hand 
side abbreviated to (mi, n2),’ 


1 pit Pre 
a(Ni, Ne) = GK) am (i iil er) (5 > na) (i oi ’ (1) 


where G(K) = G(Ki, Ko, ---, Kp) is the normalization constant so 
chosen as to make the sum of all quantities in (1) equal to 1. Explicitly, 


Bp = = pa ™ Pre 
een 2 (i (Ki - aa (3 mi) (i). ) 


m,=0 
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The function G(-) defined on the integer lattice in R” is referred to as 
the partition function.* 


2.3 System performance 


A number of interrelated system performance measures are obtained 
from the partition function. We start with N;,(K), the mean number 
of constituents of class ¢ in node 1 (‘think’), and obtain directly 
from (2), 





Nin(K) = paG(K — e:)/G(K), (3) 


where e; is our notation for the vector with the i‘* component unity 
and all other components zero. Thus, G(K — e;) is the partition 
function associated with a new population with one less constituent in 
the i" class. From (3) and Little’s theorem applied to class i and node 
1, we obtain for the throughput of constituents of class 1, 


Ai(K) = G(K — e;)/G(K). (4) 


The mean response time, i.e., time spent in the cPuU in each cycle by 
class i constituents, is obtained again from an application of Little’s 
theorem: 





té(K) = KiG(K)/G(K — e:) — pa. (5) 
Finally, the utilization of the cpu by constituents of class 1, 
u(K) 4Y ++ or m(m, Ne) = p2G(K ~ e)/G(K). (6) 
J 


The important point to note is that all the mean values given in (3) 
through (6) are simply obtained from the knowledge of the partition 
function estimated at two neighboring points of the integer lattice in 
R?. As the above quantities are all closely related, we shall henceforth 
consider only the last, {u;(K)}. 

Higher moments of N,, the random number of constituents of class 
t in node 1, may also be obtained from knowledge of the partition 
function: 


Nii(K) = (piiG(K — 2e:) + pirG(K — ei)}/G(K), (7) 
and, for i ¥ j, 
NaN y(K) = papvaG(K — e; — e;)/G(K). (8) 


Of course, the moments of {Nj2} are easily derived from moments of 
{Nix}. 


* The distribution in (1) and (2) is also the stationary distribution of other networks. 
For example, if node 2 is first-come-first-served with class independent service rate 
1/p2, then (1) and (2) with pj2 = pz is the solution. 
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2.4 Recursive solutions 


The above results explain why the problem of efficiently computing 
the partition function has excited so many researchers.’*’”*! For the 
problem at hand it is easy to arrive at the following recursion by 
established techniques: 


Nal 


PA. 


G(K) = }) pj2G(K — e,;) + J] (9) 
J=1 J= 


ext 
1 Kj! 
The boundary conditions are: G(K) = 0 for K # 0. 

Observe that the partition functions themselves can be scaled on 
account of the linearity and the fact that only ratios of partition 
function values have physical content. By the same token, implemen- 
tations of (9), for large K, typically give rise to values which are either 
very small or very large leading to rather severe problems of overflow 
and underflow. Proper scaling is only marginally helpful. 

The main problems with implementing (9) are with respect to time, 
memory required during computation, and accuracy. A straightforward 
application of (9) would require an estimated K” iterations, where K 
is the generic class size. Similarly, the storage required would be 
approximately K?~’. Now these crude estimates can clearly be im- 
proved upon by simply pruning or avoiding the computation of inter- 
mediate lattice points, but this would be at the cost of increased 
algorithmic complexity, and the extent of the accrued benefits are not 
generally known. The underflow/overflow phenomenon that affects 
the accuracy of the scheme has already been commented on; a no less 
severe problem is accumulation of round-off errors in a large number 
of iterations. 

The recursive solution in (9) is one of several that can be generated 
by recently discovered techniques. However, all solutions that we are 
aware of are recursive and share to varying degrees the three broad 
categories of limitations just discussed. 


lil. INTEGRAL REPRESENTATIONS 
3.1 Partition function 


We start with Euler’s integral 


ni = e ‘tdt. (10) 
0 
Substituting for the middle term in braces in (2) we obtain* 


* If the range of integer subscripts is not stipulated explicitly, then the range is 
understood to be [1, p]. 
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oo K K Kj-m; m; 
= —~ Fy FS Pj” (tpj2)"™" 
ace) = Jew So 3 (ge) (Ua 


m,=0 m,=0 


wus 1 sd 4 Kp Ki K; Ki-m(¢ 0) dt 
“ar J, ° 2 pj? (tpj2 


mp=0 m,=0 j 
i eee ie 
= Tz) e II (py1 + tpj2) jdt. (11) 
0 


To obtain our final form for the partition function, let N be the 
generic large parameter to be associated with a large network. Also, 
for j = 1, 2, --- p, let 


K; = B;N, (12) 
>, 4 mean think time _ ps 
” * mean service time ec PP? 
class j 
= yj;N. (13) 


The suggestion in this notation is that {8;} and {y,;} are 0(1), which is 
the situation to which our work is primarily directed. 

There is considerable latitude in selecting N. The following choice 
is certainly not essential but as it does ease some of the manipulations, 
we use it throughout the paper: 


N = ({] r®)V/2%, (14) 
An implication of this choice of N is that 
» 8; log y; = 0. (15) 


We return to (11) to observe 
Kj foo} 
G(K) = (122) | et Tu + ode 
Ki) J, 


Kj ss | 
7 (1 a) aa | a"? TI (ys + 2)" dz. (16) 
J° 0 


Finally, we have after using (14), 


Proposition 1: 


K; eo 
G(K) = (w ne) | e Nz, (17) 

J* 0 
where f(z) 4 z— > B; log(y; +z). O (18) 
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We shall see that typically there is no need to compute the term in 
braces in (17). 


3.2 Representations of mean values and some higher moments 


We have two options concerning representations of the physically 
interesting quantities in (3) through (7). Concerning ourselves with 
u;(K) given in (6), we may simply use the respective integral represen- 
tations for G(K — e;) and G(K) and obtain for i = 1, 2, --- , p, 


N | ; eNfadz 
Poo (19) 


a N i e N@Adz 
0 


where WN and f(-) are defined analogously to N and f(-) but for a 
network with one less constituent in the i‘* class. Obviously N and N 
as well as f(-) and f(-) are going to be close to each other, but the 
above option takes no further notice of this fact. 

In the contrasting option, we proceed as follows. Observe that 


i G(K + e;) 


Pi2 G(K) 2) 


ui(K + e;)7? = 


Now, from (16), 





K; 0 
= pia Pid . ae : K; 
G(K + e;) Kel (1 ) [ (r; + te I (r; + t)“idt 





K; oo 
_ __Pi2 2 Pit —Nf(z) 
= Ny vi t+ ; 9 


where to obtain the last equation we have proceeded just as we did 
from (16) to (17). The point to note is that in the above expression N 
and f(-) are identical to that used in the expression for G(K) in (17) 
and (18). Finally, from (20), for z = 1, 2,---, p 


Proposition 2: 


| ze N@dz 
1 0 


ui(K + e;)7! = yi + 7 
| e Adz 
0 


~ B+ 1/N 
Notice that the ratio of integrals, the only quantity requiring significant 
computation, is common to all classes. 


O (22) 
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The higher moments given in (7) and (8) have similar representa- 
tions which are easy to derive. We give the form for one term that 
occurs in (7) which reveals the general pattern: 

1 G(K + 2e,) 1 1 


pa  G(K) N* yi(Bi + 1/N)(Bi + 2/N) 


co co 
| ze Nladz | ze M@dz 
0 0 


ek RS rs MT (55. 


i e NAdz i e Nfedz 
0 0 


The feature to note is that the i" moment involves integrals of the 
form fo 27e/ dz, 7 =0,1, «++, i. 
3.3 Properties of the function f 
We shall later need to recall certain properties of the function f(z) 
in (18), z = 0. As a consequence of (15), 
f(0) = 0. (24) 
Note 


f(z) = 1-3 Bs/(y; + 2) fori=1 
(-1)'(i- I! VY B/Qyt+2)' fori >1, (25) 


and the following alternating sign property that holds for i = 
0: 8: 35 


1)f(z)>0, 220. (26) 
Also, for i = 2, 3, +++, 
[fl SFO], zB. (27) 


As the second derivative is positive in [0, ©), the function is always 
convex. 

Let us view the function f(z) in the interval y, < z < ©, where 
Ym 4 —minjyi (see Fig. 2). As the figure indicates, the derivative of the 
function tends to o at both ends of the interval. Coupled with the 
convexity and the other already established facts, it follows that the 
function has a unique stationary point, a minimum, in (ym, ©). As in 
the figure, we let Z denote this unique stationary point. 

This stationary point may be obtained as the unique real solution in 
(ym, ©) to the equation 


B; 


a ytz 





= 1. 


Me 


Thus, the largest real root of a polynomial of order p gives Z. 
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F(z) f(z) 





Fig. 2—Sketches of the function f(z). 


It will be important for us to distinguish the cases Z < 0 and Zz = 0. 
The slope of the function f(z) at the origin effectively indicates which 
case holds: the first if the slope is positive and the second otherwise. 
Thus, a key parameter of the system is 


a & —f%(0) =¥ Bi/y-1 
=) Ki/r; — 1. (28) 
From the previous discussion we thus have that 
Min f(z) = f(0) if as0 
= f(z) if a>0. (29) 


Equation (29) summarizes the background on the stationary point 
as needed for most of the paper. For example, Z is needed only if a > 
0. However, Section VIII is exceptional in that, while considering small 
possibly negative a, the corresponding stationary point Z is required to 
be known. Note that as a — 0, 2 ~ a/Y B:/y?. 

The parameter a, a > —1, is an indicator of the traffic intensity, 
with increasing a corresponding to higher traffic intensities. Our re- 
sults, theoretical and numerical, show that a = 0 corresponds to heavy 
usage corresponding to close to 100 percent utilization of the cpu. 
‘Normal’ usage in large networks will certainly require a < 0 and in all 
likelihood a will not be close to 0. For this reason, the most compre- 
hensive results given here are for the case a < 0. 


IV. ASYMPTOTIC EXPANSIONS FOR NORMAL TRAFFIC 
Throughout Section IV we shall consider a < 0. 
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4.1 Laplace’s method 


We shall first apply Laplace’s method” to obtain asymptotic 
expansions (see Appendix A for notation) for the integral in the 
representation of the partition function in (17), and subsequently to 
the other integrals in (22) and (23). Laplace’s method observes that 
when NV is large, the minimum of Nf(z) at z = 0 is very sharp and that 
the dominant contribution to the integral comes from the neighbor- 
hood of 0. 

In the integral 


I= [ eM az. (30) 


0 


let us change variables, 


u& f(z), (31) 


= 7 —Nu dz 
[= [ e (2) du. (32) 


To obtain dz/du, let us begin with the power series convergent in the 
neighborhood of 0, 


so that 


u= f(z) = x fiz, (33) 
j2 


where f; = f”(0)/j!. Standard procedures allow the series to be 
reversed, i.e., give Zz aS a power series in u. Appendix A elaborates on 
this procedure and explicitly gives the leading terms. Let the series 
obtained by this procedure be 


zAg(u) = ¥ gy’, (34) 
J21 
so that 
dz j 
aan » aju’, (35) 


where a; = (j + 1)gj+1. Substitution of this series in (32) yields 


I~) a | e N“y/du 
J2=0 0 

and thus 

Proposition 3: 


i 
faye 2 Noe: ‘Oo (36) 


j=0 N’* 
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The leading three terms of the series are obtained from 
a=1/fy  a=—2fe/fis a2 =3(2f3— fifs)/fi. — (37) 


The proof that (36) is an asymptotic expansion is available from 
various sources.’*””*? Indeed an explicit proof is available from the 
bounds that we develop in the following section in the course of an 
error analysis. Nonetheless, we sketch a proof based on Watson’s 
lemma™ applied to the integral in (32)—the lemma is anyhow used 
later. Watson’s lemma considers the integral {¢ e%“h(u)du in which 

(t) h(u) has a convergent power series expansion in the neighbor- 
hood of the origin, and 

(ti) there exist constants ci, c2 such that | A(t) | <cie*’ for ¢= 0, and 
asserts that an asymptotic expansion for the integral is obtained by 
replacing h(u) by its power series and integrating term by term. Thus, 
series reversion to obtain dz/du and a subsequent application of 
Watson’s lemma gives the asymptotic expansion in (36). 

The reader will note for future reference that the asymptotic expan- 
sion in (36) may also be written as 


I~ ¥ g(0)/N?, (38) 
j=l 


where g(-) = f~*(-), as follows from (33) and (34). 
Asymptotic expansions of integrals of the form I = {[¢ z*e~*/dz 
follow with only slight modifications. Thus, in lieu of (32) we now have 


r= | om (eS) du (39) 


0 


In principle, it is straightforward to obtain a power series for (z“dz/ 
du), which is convergent in the neighborhood of 0 from the power 
series in (34). Using the following as the defining relation* for the 
sequence {a\*’}, 7 = 0, 1, 2, ---, 


k 
zt (x ew’) (x (i+ Deve’) = Yaj’u’**, (40) 


du jz jz0 j=0 


the asymptotic expansion for the integral is obtained after term-by- 
term integration, giving 


Proposition 4: 
I® = | ze MOdz~ Y (j + k)lay?/NT (41) 
0 j-0 
as N>-o, O 
* In this notation, the sequence {a;} in (35) and (36) is {a;”}. 
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Let us consider in greater detail the expansion of the integral J‘ = 
fe ze—*’ dz, which will be needed if (22) is used to compute the mean 
values. We have derived the following recursive formula which effi- 
ciently generates {a 1 from the coefficients {a;”} that are needed 
anyhow for the expansion of J: for all 7 = 0, 


-5 AS nai Ah ik. (42) 


In particular, the leading three terms of the expansion of J" are 
obtained from, 


Y= 1/fi; af? =—-3f/fi; af? = 2(5f3 — 2Af)/fe. 


We have derived, but omit to give, a more general recursive for- 
mula—of which (42) is a special case—for generating {a{**"} from 


{a;"}. 


4.2 Asymptotic expansions for the utilizations 


As discussed so far, both of the two options stated in Section 3.2 for 
representing the mean values [see (19) and (22)] require the develop- 
ment of asymptotic expansions of two separate integrals and the 
subsequent computation of the ratio. Here we observe that a single 
asymptotic expansion exists for the mean values. The un- 
derlying reason is that the asymptotic sequence {N~¥}, j = 
0, 1, --- , form a multiplicative asymptotic sequence.'*” 

In particular, if as in (36) and (41), 


2 (0) (0) (0) 
= ao ay A> 
| e M@dz ~—__+—, + M+. 
0 


‘N2 N 
and 
oo (1) (1) (1) 
| ze N@qz ~ Go Bee qe vee 
5 N N N 
then 


| ze N@dz 
0 bj 
=a 
e Nledz 
0 


where the sequence {0;} is obtained by formal substitution. 
The following gives the leading terms for the utilizations derived by 
the above procedure. 
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Proposition 5: 
Fori=1,2,---,p,asN—> %, 
Ai Aso. As 


{ui(K + e:)}" (Bi + 1/N) ~ y+ + 5+ Ne 


NN? ) 


where 
Ai=—l/a; A2=4f/a*; As = —40f3/a° — 18f/a*. O 


In accordance with an earlier observation, all terms of the asymptotic 
series other than the first are independent of 1. 

Proposition 5 contains a justification for treating the parameter a 
[see (28)] as an indicator of traffic intensity. Using only the dominant 
term, 


{u;(K + e:)}~* ~ yi/Bi (44) 
and 
utilization of cpu = )' u;(K) ~ 1+ a. (45) 


A necessary caveat is that the above, as indeed all results in this 
section, has been derived for the assumption a < 0. 

Since the utilization as given by (45) can come close to unity even 
with a < 0, (45) justifies another earlier statement that for large 
networks normal usage will not extend beyond the range a < 0. 


V. ERROR ANALYSIS AND PERFORMANCE BOUNDS FOR NORMAL 

TRAFFIC 

Maintaining the restriction a < 0 placed in the preceding section, we 
supplement in two directions the results obtained so far. First, we 
obtain essential results on the error incurred in truncating the expan- 
sions. These results containing information extending beyond what is 
required as proof of asymptotic expansions are needed for very prac- 
tical reasons, such as to know how many terms to use and, more 
importantly, to help define the regime of applicability. In the second 
part of the section, certain rather special properties of the functions in 
the representations are used to derive analytical bounds on the net- 
work performance measures. 


5.1 Completely monotonic functions 


The following result on the function g(-) = f~'(-) is key to much of 
the error analysis: 


Proposition 6: 
(-1)/’g(u) <0 for u=0, j=1,2,3,--- O (46) 
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An inductive proof is given in Appendix B. By virtue of this result, 
g'(-) isa completely monotonic, or alternating, function (see Ref. 23). 
The importance of this property stems from the role of g“(u) = dz/ 
du in the integral representation (32). 


5.2 Error bounds 


In connection with the expansion of the integral J in Proposition 3, 
let R,, denote the error that accrues if only the leading m terms are 
used, 1.e., 


m-1 m 
Rn=I- ¥ jlaO/N*=1-— Y g%(0)/N4. (47) 


j=0 j=l 


Now by the mean value theorem,” for each u there is a £(u) in [0, uw] 
such that 


gu) = x g(0)ul"/(7 — I+ go" (G)u"/m! (48) 
jc 
On substitution in (82), 
I= 5 g/(0)/N’ + Rm, (49) 
= 
where 
Rn = = ; e Neg ™ Dr e(u)}u™du. (50) 


A simple corollary to (50) and Proposition 6 is 


Proposition 7: 


Rn>O if mis even, 
<0 if misodd. UO 


This, of course, means that the terms in the asymptotic expansion 
alternate in sign and that the partial sums of the asymptotic expansion 
alternately over- and underestimate the true value of the integral in 
the following manner. 


Proposition 8: For m even, 


7 i=] 


i= 


m id m+1 
y g(0)/N'< | edz < F gO) /NL O 
: 0 


The above is quite useful since in most situations the designer would 
much rather overestimate than underestimate a measure such as CPU 
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utilization. In this context, both the upper and lower bounds are 
required since ratios of integrals occur in the measures. 

An implication of Proposition 6 is that |g*(é)| < |g*”(0)|, 
£> 0, which together with (50) gives 


Proposition 9: 
Rn < g*t?(0)/N”™*' if mis even 


> g™V(0)/N™*' if misodd. O 


The above propositions thus state that the error is numerically less 
than the first neglected term of the series, and has the same sign. In 
particular, we have an explicit proof that (36) constitutes an asymptotic 
expansion. More generally, the above results show that, on account of 
the specialty of the integral, the main results that we require from an 
error analysis are already present in the expansion. 

It is useful to examine in detail g“ (0) and thence the bound on Rs: 


|Rs| S (bs|a|? + 10a] b2b3 + 1503) /(la|7N%), (51) 


where }; = (i — 1)! ¥| B;/yi. (A look at the proof in Appendix B will 
convince the reader of the presence of |a|?”*'N”*? in the denominator 
of the bound for |R,|.) The bound does make the suggestion that in 
cases where a is so small [i.e., utilization is very large, see (45)] that 
aN is itself small, then |Rs3| is large. More generally, in the case of 
small «”N, the number of terms in the series requiring computation to 
meet specifications on the accuracy may be large. Later we return to 
consider this case further. 


5.3 Bounds on mean values 


The following bounds which supplement the computational proce- 
dures are presented to serve as aids in design and synthesis. For z = 1, 


2, eee) D, 
~ Proposition 10: 


| (we +e)}7 are (B; + 1/N) 


Bit 1/N 
oie, V1 + 8f/a°N (52) 
7 2|a|N 
eB) - ee C 53 
—2hN\ 1+ V1 + 16fa/(7a?N) i 
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Recall that 2f = f(0) = ¥: B;/y? while a (here a < 0) and N are as 
given in (28) and (14). 

To prove Proposition 10, we may use the representation in (22) for 
u; in which case it remains to bound from above and below the pair of 
integrals appearing there. This is done in Appendix C by making use 
of the sign properties of the higher-order derivatives of the function 
f(-). | 

Notice that as N — o, the upper bound in Proposition 10 on 
{u;(K + e;)}~' approaches {y; — 1/(aN)}/(B; + 1/N), the sum of the 
leading two terms of the series for u; in (43). Also, as N — o, the 
expressions in (52) and (53) both approach 0. 


Vi. EXPANSIONS FOR THE CASE a = 0 


At the end of Section 5.2, we commented that the expansions given 
earlier may require a large number of terms in regions where « = 0. 
For this case, we here generate a somewhat different and more efficient 
series. We shall, however, be quite brief in our exposition because the 
uniform asymptotic expansions to be derived in Section VIII also allow 
appropriate expansions to be obtained. The ad hoc but direct treatment 
here is supplementary. 

We need notation that is specific to this section. Let 


K; = iN + dWN 7 
peGNECIN oor (54) 


Each of the variables {d;} and c may be either positive or nega- 
tive. However, as our interest in this section is in a ~ 0 and as a ~ 
¥} bi/ai — 1 as N > ©, we require that 


> bi/ ai = 1, (55) 


In a computational procedure the above restriction poses no particular 
problem. 
Also, we shall mainly consider the integral 


TA (]] rj *) | e ‘TT (n+ t) “dt. (56) 


Comparison with (16) shows that the integral is related to the partition 
function thus: 


K; 
G(K) = (1 a (57) 


As previously, the computation of the quantity in parentheses is not 
required to obtain the mean values. Our main result is 
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Proposition 11: As N > », 
fe DG NNT, (58) 
j=0 


where the sequence {c;} is given below. UO 


The proof of the proposition is in Appendix D. Here we comment 
on some features of the sequence: 


C= | e~40°/2-Br ET (py) dy, (59) 
0 


where A = ¥\ b;/a? > 0 and B = c — ¥. d,/a; and H;(») is a polynomial 
of degree 37 in v with coefficients that are fairly straightforward to 
obtain. The key point is that A, B as well as the coefficients of the 
polynomials are all 0(1), i.e., N does not enter into their definitions. 
Results given in Section VIII indicate how the coefficients c; given in 
(59) may be effectively computed. 

We give below the leading three polynomials in the sequence 
{H;(v)}: 


Ad(v) = 1, 
H,(v) = ¢ — (¥, dj/a?)v?/2 + (> bi/a})v?/8, 
H2(v) = —c( > d;/a?)v?/2 + (¥. d;/a? + c ¥ bi/a?)v3/3 
+ [(¥ d;/ai)? — 2 ¥ b;/af]v1/8 
— (¥ di/ai)(¥ bi/a?)v°/6 + (¥ bi/a?)?v°/18. (60) 


Vil. ASYMPTOTIC EXPANSIONS FOR HEAVY TRAFFIC CONDITIONS 


Here we obtain asymptotic expansions for the basic integral in (17) 
and (18) for the case a > 0. For reasons similar to those discussed 
earlier for the case a < 0, the expansions to use for a ~ 0 are in 
Sections VI and VIII. Hence, the expansions given below are for 
exceptionally heavy traffic conditions, where a is not only positive but 
also not close to 0. 


7.1 Laplace’s method 


The key difference from the treatment in Section 4.1 is the presence 
of the singularity at Z (see Fig. 2), which will be assumed to be known. 
For large N the dominant contribution to the integral 


r= | e N@qz (61) 
0 
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comes from the neighborhood of Z. A Taylor series expansion around 
Z gives 


f(z) —fo= x f(z - 2), (62) 
where 
f=f%A/, F=0,1,2, +++. (63) 


In particular, for 7 = 2, 
f= (-1) b Bil (yi + av Ji (64) 


We make the following specific decomposition of J, which is conven- 
ient here and even more so in the error analysis to follow with ze as in 
Fig. 2, 

2 7 22 ‘ 
eNhy = | e NUl@)-h} dz + | e Nf@-h} do 


0 2 


+ | e NUlA-h gy, (65) 


2 
Consider the terms in turn starting with the middle term. If we let 
udf(z)-f, 224 (66) 


and use the series in (62) for the right-hand side, then we may reverse 
the series, as discussed in Appendix A, to obtain 


foo} 


z-ZAg(u) = ¥ gw” (67) 


j=l 


with the coefficient g; depending only on the coefficients fy, fs, --+, 
f;. Now, 


20 . —fo 
| e Nh) ds = i ea Ou)du 
Ps 0 


z 


—fo [ove) 
= | ee SY aul du, (68) 
0 


j=0 


where a; = (j + 1)gj+:/2. The individual integrals in the sum will be 
recognized to be incomplete gamma functions. 

On returning to (65) and the first term in the right-hand side, we 
find by an identical argument that 


618 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1981 


z 2 fo oo 
| oe NU-A) da = ey (—1)/aju du. (69) 


0 0 j=0 
The two integrals in (68) and (69) may conveniently be combined to 


give 


—fo fo) ad 
eviy = | eX" SY Qaojul du +| e NU@-fh ga, (70) 
0 J=0 z. 


2 


At this stage, the following two approximations are made, with their 
effects bounded in Section 7.2 in the course of the error analysis: the 
integration interval in the first term is extended to [0, ©) and the 
second term is ignored. Nonetheless, the error analysis shows that 


I~ er | eX” YS Qaoju’ "du, (71) 
0 g=e 


giving Proposition 12. 
Proposition 12: As N— o, 


I= | e M@dqdz ~ e-Nh YX 2B + Y)aa/N. O (72) 
0 j=0 


Recall that (4) = Va and for7 =1,2,---, 
T'(j + 4) = Vz [fin G— ¥). 
We give the leading three coefficients: 
ao = (4) /f2”, 
a2 = (Ye) (15/3 — 12f2 fa)/F2”, 
as = (5456) (—64f3fs + 224f3fofs + 112f3f% 
— 504ff3/, + 231F4) /f2P”. (73) 
The procedure for obtaining the asymptotic expansion for the 
integral 
IY = | ze W@qz (74) 
0 


is similar. Notice 


d® — ZI) = em | (z — Ze NU@-h) gz. (75) 
0 
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we find it convenient to expand the integral on the right-hand side. 
The expression analogous to (71) is 


foo} 


(I — 21) ~ e-Nh | eN Y) 2a3)-.ul?du, (76) 


0 Jai 


where, for 7 = 0, 1, 2, ---, 
1 J+1 
ay? == Y) Rgj-n+28k, (77) 
2 k=1 


which is to be compared with the expression for a; following (68). The 
following is a useful formula for efficiently generating the sequence 
{a9} from {a;}, which is needed anyhow for computing I: 


J 
a? =2 Y arra/(k+1), f= 0,12, ++. (78) 


Recognizing the gamma functions in (76) gives Proposition 13. 


Proposition 13: 
, 2 1 : 
I — ZI ~ e7N* ¥ ar (i * 5) a2ja/NP? (79) 
j=l 
as N-> », 0) 


The asymptotic expansions in Propositions 12 and 13 may also be 
combined, as discussed earlier in Section 4.2 to yield an asymptotic 
expansion for the mean values. In particular, we obtain Proposi- 
tion 14. 


Proposition 14: As N—> ©, 
(u(K +00) (8+ Y/N) ~ + 44S, 
N N 
where 
Ai = -3fs/4f3, 
Aa = (6 fafa fs — 15 f3fs/8 — 135 f3)/f8. 0 (80) 
7.2 Error analysis 


The analysis to be presented supplements the result in Proposition 
12 and the error estimates to be given provide guidelines for the use of 
the expansions. The broad outline of the analysis have been suggested 
in Ref. 23. 

As in Section 5.2, let R, denote the error incurred when only n 
leading terms of the series in Proposition 12 is used, i.e., 
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. n-il 1 ; 
Ri =I-e% 7 ( jt 5] Qar;/NI*Y?, (81) 


j=0 


For I we will use the expression given in (70) and decompose the error 
R, thus 


Rn = —é€na(N) + €na(N) + €3(N), (82) 
where 
m 2 n-1 
€ni(N) = ene | e Nu ( YX 2a") du, (83) 
=f, j=0 . 
. —fh foe) ; 
€n2(N) = em | eo > 2a!) du, (84) 
0 =n 
€;(N) = | e Nf)dz, (85) 


2 


Thus, the three terms on the right-hand side of (82) respectively 
denote components arising from the extension of the integration inter- 
val from [0, —fo| to [0, ) in (70) and (71), the use of only n leading 
terms from the infinite series in (71), and the neglect of the second 
term in the right-hand side of (70). Each component is now bounded. 

To bound ¢,,:(N)), we make use of known bounds on the incomplete 
gamma function: 


| e Neysi-V2qy =[ (i 4 4 = wh) [wee 


-hy 
eNh(N| fol)? 


ae 
weal Nf — max (i- > 0) 


Thus, 


ae | 
for n|f| > max (i-}0). 


2 = 
€n,(N) | = —————_ Q; he 86 
| én (V) | Nfl ae by (ee! ll (86) 
where 6, = max(n — 3/2, 0). 
In bounding €,,2(N), we will postulate the existence of a finite valued 
o, with the property that 


le 2) 
> 2a2 ju? /? 


jan 


< |2a2,|u"-/e",0<u<| fol. (87) 








This approach fails when o, is infinite but the characteristics of the 
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integral and the specific decomposition (65) that has been employed 
preclude this possibility. Let 


On = max y,(w), (88) 
lul<l fol 


where 

iv o} 

Y, 2aeju7-"”? 
j=n 


89 
2aontt oo) 


n—-1/2 








won 
u 


Small u is where (see Ref. 23) the danger of unbounded o, is usually 
most manifest. However, as 


2 

Q2n+2 a a 

SO epi ope gece) (90) 
Q2n Aan Aan 


when u — 0°, no problem arises here. 
Using (87) in the defining expression for €,2(N) in (84) 


I fol 
| €n,o(V) | = e Nho| Qaan| [ eo N-onu 2-1/2, 
0 


= e Nhl Qa, |T(n + %)/(N — on)", (91) 


The last term to be considered from (82) is e3(N). We use the 
following property. 


f(z) = f'(z)(z— 22), 22 22, (92) 
which yields 


les(N)| Ss (93) 


1 
Nf’ (z2)’ 
a small quantity compared to the right-hand side of (91). This con- 
cludes the process of bounding the components of the error term R,. 
A corollary is the proof to Proposition 12. 

The bound in (91) is the largest component in the error bound. In 
examining (91), we observe that the condition in which the bound is 
large is when a2,/N” is large. Now the expression for a2, contains in 
the denominator a term /3”*"””, as (73) attests. Thus, whenf3NV is small, 
we expect the asymptotic expansions in Section 7.1 to be inefficient. 
We return to this case later in Section 8.2, where this as well as the 
similar difficulty encountered in Section 5.2—where a was negative 
and small—is treated in a unified manner. 


Vill. UNIFORM EXPANSIONS 


This section has two objectives. The first is to show that there is a 
framework that unifies the expansions in Sections IV, VI, and VII. 
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This consists of showing that the integrals of interest may each be 
given by a common expansion valid uniformly for the entire range of 
values of the system parameters. These expansions turn out to be in 
parabolic cylinder (or Weber) functions.’*”°”* The advantage derived 
is that these classical special functions have extensively documented 
expansions for the entire range of parameter values.”*”* Indeed, using 
these expansions we sketch in Sections 8.3 and 8.4, at the cost of some 
duplication, derivations of the expansions obtained earlier in Sections 
IV and VII for a < 0 and a > 0, respectively. The second objective is 
to derive a computationally efficient expansion for the case a ~ 0, Le., 
where the stationary point Z is very close to the boundary of the 
integration interval. The error analysis in Sections 5.2 and 7.2 has 
shown the need for a separate treatment. The expansion that is 
obtained for this case in Section 8.2 is obtained from an appropriate 
expansion of the parabolic cylinder functions. 


8.1 Uniform expansions in parabolic cylinder functions 


Consider the integral 


I= ‘ e Nl2qz (94) 
0 


without restrictions on the parameter a. Following Friedman,” con- 
sider a change of variables from z to v given by 


v? — 2avu = f(z), (95) 


where a is a parameter of the transformation to be fixed later. The 
objective of the transformation is that the component of the integrand 


‘in braces below 
I= | as, (=) dv (96) 
, dv 


satisfy the dual requirements of boundedness and a convergent power 
series, as required for an application of Watson’s lemma (see Section 
4.1 following Proposition 3). [The reader may verify that the simpler 
transformation v = f(z) violates the boundedness requirement when- 
ever a > 0 and z = Z since f (2) = 0.] For the transformation in (95), 


dz _ 2(v-a) 
du f(z) : 


This suggests the selection of the parameter a to be such that v = a 
when z = Z, with the accompanying indeterminacy and the possibility 
of boundedness of dz/dv. This key clue does indeed give a unique map 
of the form in (95) with the desired properties, as summarized below. 


(97) 
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Proposition 15: For z = 0, let 


u(z) = a + sgn(z — 2)V f(z) — F(Z), (98) 
where the constant a depends on all the system parameters: 
a = (sgn a)Vv—f(2). (99) 


The transformation is monotonic, increasing and maps [0, ©) to 
[0, 0). Also du/dz is continuous and uniformly bounded. U0 


The transformation is used to derive for dz/du a convergent power 
series Yo C,v’ in a neighborhood of the origin. This is achieved in three 
steps: 

(t) Use (98) to obtain 


v(z) = Ayz + Aoz? + Agzi t+ ---. (100) 
(tt) Reverse the series (see Appendix A) to obtain 
2(v) = By + Bo? + Bgv? + ---. (101) 
(iii) Differentiate term by term to obtain _ 
d 
<= Cot Cw + Co? +o, (102) 


where C; = (7 + 1)Bjs1. 

The reader may verify that the leading terms of the sequence {C;} 
thus obtained are as follows [recall from (33) the definition f, = 
f' (0/711: 

Co = 2a/ Qa, 


2 
a= =| (*) p-1| (103) 
Qa Qa 


The above tacitly assumes that 2 and hence a have been evaluated. 
The power series expansion for dz/du may now be substituted in 
(96) to yield 


I= SG | e Nw? 2a) nid, (104) 


J=0 0 


The integrals appearing above are simply related to the parabolic 
cylinder functions U(., -); thus: for 7 = 0, 1, 2, ---, 


. —N(v2—2av) Jd —_ eNa?/2j : 1 J2N 
; e U v= Onan U DG 2N |. (105) 


Expansions for related integrals such as 
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I) = i ze Nladz (106) 
0 


are only slightly different. Here the term dz/dvu is replaced by zdz/du, 
which has the power series expansion 


zdz/dv = ¥. C#v!, (107) 
j=l 
where, see (101), 
j-1 
CP = ¥ B)-2Cr. (108) 
k=0 


Specifically, 
ow = (22) 
gece 
2 
or-2(2)((a)e-n} 
ala a 


The following summarizes the expansions in parabolic cylinder 
functions of the two integrals of main interest, with the expansion 
valid for all a. 


Proposition 16: 


= = = ? ; 
| e N@dz= ¥ a| e Nv 2a jdy 
0 0 


j=0 
— pNa? . IG; 7 1 
or Dany ere u(j ga ag avi), (110) 
I? = | ze Nfadz 
0 
— ,Na?/2 is TCP . dl JAN 
= eM" Youn U(s+5.-av2N}, (141) 


where the sequences {C;} and {C9} are respectively as obtained by 
the procedures in (100) through (102) and (108). Specifically the 
leading terms are as given in (103) and (109). O 


We should add that the above expansions are not strictly asymptotic 
expansions since the parabolic cylinder functions do not satisfy the 
requirements of asymptotic sequences for certain ranges of the param- 
eter a?N.* The interested reader will find in Ref. 27 a description of 
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the process for obtaining uniform asymptotic expansions of the inte- 
grals. However, we have not found it necessary to undertake the 
additional effort required to obtain the coefficients of the uniform 
asymptotic expansions. This is because for aN small, the case treated 
below in Section 8.2 and of main interest, the functions U(j + %, — 
av2N), j = 0, have all the desirable properties that are required of 
asymptotic sequences. 

A noteworthy property of the functions U(-, -) that can be important 
in computations is that it satisfy the recursion” 


xU(j + %, x) = UJ — %, x) — (7 + IU + %, x). (112) 


8.2 Expansions for the case of a?N small 


To motivate the results to be given here, observe that the stationary 
point of the curve of f(z) (see Fig. 2) is close to 0 when a is small 
(utilization of cpu is high), since 


Z~a/(2fr) ~ as z—>0. (113) 
On enquiring how the parameter a behaves for small Z and a, we find 


from (99) that 


a =——.,. + O(a’) as a— 0. (114) 


Al 


Thus, the case of small a, which we know from Section 5.2 requires 
special treatment, corresponds to small a. 

Small a is also implied by small f. and is therefore also of interest on 
account of the discussion in Section 7.2. This follows from 


a = 2( fo)? + O(2°”) as z—>0. (115) 
For small a?N, it is known” that for j = 0, 


-Nw?-2av) nig = Nara ul; ee — avaN 
: c i oN) FD? J 2’ 
hid GV) SAT) (A) 26) 4 
~ NOW? [uo + (aVN) ui?’ + (avN) "p27" + +++], (116) 


where 
noe J ie ia 
i! T(1 + j/2) 
y! 9(+1)/2-j-1 
aT + D2) 


In these relations, notice first in (116) the desirable presence of the 
powers of N in the denominator. Secondly, in connection with the 


(7+ 1)(7 + 3) --- (JF +7t—-D), leven 


(7+ 2)(7 +4) --- G+i- 1), iodd. (117) 


626 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1981 


sequence {!”}, observe that p{Z}/ul) = 2(7 + i + 1)/{(i + 1)(i + 2)}. 
Thus, for fixed 7, the sequence converges rapidly to 0 with increasing 
i. These observations state that when using the expression (116) in the 
expansions of Proposition 16, first, it is necessary to compute only 
up to small values of the index j and, secondly, with aVN small, the 
computation of the bracketed quantity in (116) also needs very few 
terms. 

The following summarizes the important computational procedure 
described above. 


Proposition 17: For small one 


T= vn z eee b (aVN)'n | (118) 
(1) 
[= Pe Dy —— | x (avy? | (119) 
Jj21 = 


where pi!” is in (117), and a, {C,}, {C$} appear in Proposition 16. 0 


8.3 Expansions for normal traffic - 


For normal traffic, a < 0 and consequently a < 0. If in addition, 
a <0 or, specifically, a?N >> 7”, then™ 


| e Nw?-2ar) Jy 
0 
1 4» (+2) (7 +4)! 
ee ESSE (a Ents oAmnash on Span eeu, rs) NEB ?, §) 
(—2aNn)™ ( 4(@=N)  32(a°N)? a) 
It can be shown after some manipulations, which we omit, that this 
expansion when substituted in Proposition 16 is identical to the main 
result of Section IV, namely, the expansion in Proposition 3. The 
bridging relation is 
(v—a)?=uta’, (121) 
where u is the integration variable in Section IV [see (31)] and vu is the 
similar variable in the uniform expansion [see (98) ]. 
8.4 Expansions for heavy traffic 


For heavy traffic, a > 0, and therefore a > 0. When a?N > 1 as well, 
then” 


| e N(v2-2av) dy a one NT; 1 AS ian 1) 
: JN 4a°N 


(J-1(7-2)(7-3 
aca Bec -) (122) 
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Notice the departure from the expressions in (116) and (120) in the 
absence of N’ in the denominator. 

It may again be established, although not in a simple manner, that 
the main result of Section VII, the expansion in Proposition 12, is also 
obtained by substituting the above expansion in Proposition 16. The 
bridging relation is 


(v-a)’ =u, (123) 


where the variable u is as in (66). 


IX. COMPUTATIONAL NOTES 


The asymptotic expansions of integral representations of various 
quantities in the closed queueing networks discussed above, have been 
implemented as a user-oriented interactive package on Digital Equip- 
ment Corporation’s VAX 11/780 operating under programmers work 

_ bench UNIX* system version 3. The package written in C-language 
has about 400 C-language statements and occupies about 60 Kbytes of 
storage. 

The number of classes and constituentst that the package can 
accommodate is so large that in effect no restriction is placed on these 
parameters. The other main features of the package are enumerated 
below. 

(t) The package is user oriented and easy to use. The user is 
prompted for relevant problem data. As this is supplied, validation and 
feasibility checks are made and the user informed of any errors. 

(ti) The output of the package includes all relevant statistics on 
each class (response time, utilization, etc.) including the percentage 
error incurred in the expansion. As an option to the last named, the 
user can display the various terms used in the expansion. 

(tit) The package is partly adaptive in that it automatically detects 
the divergence of the asymptotic expansion and truncates the series at 
the point of divergence. 

(iv) Numerical stability is enforced by proper choice of N. 

Numerous computational experiments were performed to compare 
the efficacy of our package, called ASYM, to the current version of a 
popular commercially available package CADS. CADS is marketed by 
Information Research Associates and a version of it runs on a VAX 
11/780 operating under UNIX time sharing system. The test problems 
run on both these packages are real-world problems encountered in 
performance analysis of a Bell System project. 


* Trademark of Bell Laboratories. 
{In the computer science applications, the class sizes in closed networks {K;} are 
referred to as degrees of multiprogramming. 
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The results of the experiments indicate that CADS is unable to 
solve our moderately sized, closed-queuing-network problems. Two 
features that accompany the breakdowns are noteworthy. One is 
numerical instability which manifests itself as overflows or underflows. 
Recent experience is that rescaling the service rates, an often-men- 
tioned device to combat the problem, does not help in substantially 
increasing the problem size. This accuracy problem is of course relieved 
with the use of more powerful floating-point machines like the CDC 
6600 (manufactured by Control Data Corporation) or IBM 3033. The 
other feature is the built-in limitations of the package. For instance, 
the current version of CADS limits the degree of multiprogramming 
for any one class to 100. There is also a limit of three classes. 

The above discussion is not intended to disparage the usefulness and 
success of CADS. CADS is extremely powerful in solving small 
queueing-network problems. Packages implementing recursions, as 
CADS does, and implementing expansions of integrals complement 
each other and, when integrated, will provide a powerful general-. 
purpose package. 

The computational experiments (see Tables I to IV below) yielded 
the interesting fact that the asymptotic expansions are quite effective 
(i.e., yield a small percentage of errors) even for small problems, 
provided the right expansion is used. This in conjunction with the 
linear growth in computation time with increased number of classes 


Table I[—Results for test problem 


No. 1* 
Degree of fe gt Utilization of 
Multipro- af srecrals a oe CPU Given 
gramming y by ASYM 
10 0.0417 0.0414 
20 0.083 0.0829 
30 0.124 0.124 
40 0.166 0.165 
50 0.207 0.207 
60 0.249 0.248 
70 0.290 0.289 
80 Breakdown 0.331 
*90 Breakdown 0.372 
100 Breakdown 0.413 
110 Breakdown 0.618 
150 Breakdown 0.618 
200 Breakdown 0.819 


* Problem specification: 
No. of classes = 1 
Think time = 240 seconds 
Processing time = 1 second 
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Table Il—Results for test problem Table IIl(a)—Problem specification for 





No. 2* test problem 3 
Total Utilization Total Utilization ‘ j Degree of 
eee Oa, of CPU Given by of CPU Given by Service Rate Service. “Multipro- 
CADS ASYM Class Server for CPU for oror This” 
10/10 0.118 0.119 This Class This Class Class 
20/20 0.239 0.23 
1 0.0001 6.00 500 
30/30 0.358 0.35 2 0.0005 2.00 500 
40/40 0.476 0.48 
3 0.0006 4.00 500 
50/50 0.593 0.60 4 0.0002 0.33 500 
60/60 Breakdown 0.70 
5 0.0002 0.60 500 
80/60 Breakdown 0.72 6 0.00005 0.20 500 
90/60 Breakdown 0.97 
7 0.00005 1.00 500 
100/50 Breakdown 0.69 
110/50 Breakdown 0.71 
140/50 Breakdown 0.79 Table IIl(b)—Results of test problem 3 
200/10 Breakdown 0.54 output by ASYM 
170/40 Breakdown 0.75 Response : 
i ea. ate % Error in 
* Problem specification: Class ( a s) Utilization Utilization 
No. of classes = 2 
Think time for class 1 = 450 seconds 1 0.95 0.01 0.0025 
Think time for class 2 = 150 seconds 2 2.85 0.12 0.0025 
Processing time for class 1 = 1 second 
Processing time for class 2 = 1.5 seconds 3 1.42 0.08 0.0025 
4 17.21 0.3 0.0025 
5 9.48 0.17 0.0025 
6 28.46 0.13 0.0025 
7 5.70 0.02 0.0025 


Table IV(a)—Problem specification for 
test problem 4 


Class 


Od Oo PO Nr 


11 
12 


13 
14 


15 
16 


17 


Service Rate 
of Infinite 
Server for 
This Class 


0.0033 
0.033 


0.0033 
0.033 


0.033 
0.033 


0.033 
0.00033 


0.00055 
0.0033 


0.00033 
0.00055 


0.0003 
0.033 


0.00033 
0.00055 


0.0003 


Service 

Rate of 

CPU for 
This Class 


to ) 
SO PBR WS 
bb NN AD ra) 


rire ro Se fe 


Degree of 
Multipro- 
gramming 
for This 
Class 


a aan aan oo om oOo oN aa ao 


Table IV(b)—Results of test problem 4 
output by ASYM 


Class 


oy Oo BO NF 


li 
12 


13 
14 


15 
16 


17 


Response 
Time 
(seconds) 

0.09 
0.834 


0.43 
0.42 


0.42 
0.28 


0.09 
2.85 


2.9 
2.82 


8.530 
8.523 


8.540 
1.63 


1.71 
1.72 


1.71 
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Utilization 
0.008 
0.080 


0.004 
0.046 


0.040 
0.027 


0.008 
0.003 


0.005 
0.027 


0.008 
0.013 


0.007 
0.156 


0.002 
0.003 


0.001 





% Error in 
Utilization 


0.05 
0.045 


0.05 
0.048 


0.049 
0.049 


0.05 
0.06 


0.05 
0.05 


0.05 
0.05 


0.05 
0.04 


0.05 
0.05 


0.05 
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makes the use of the asymptotic expansions attractive even in cases 
where the recursive implementation does not break down. 

Tables I and II display the results of problems run both on CADS 
and ASYM. Tables III and IV show the results on two large problems 
that were not admitted by CADS, but were solved with good accuracy 
by ASYM. 

It may be observed from Table I that the results from ASYM and 
CADS agree rather well, even for small degrees of multiprogramming. 
In the cases solved by both CADS and ASYM, a was small and N 
(about 240) large. Observe also that CADS was unable to solve cases 
with K;, larger than 70, even though higher values of K, correspond to 
quite low usage of the cpu and are quite interesting. 

Table II shows that CADS also broke down on a relatively small 
problem with 2 classes and 60 customers in each class. 

Tables III and IV display the output given by ASYM for the large 
problems. Table III corresponds to a problem where the total degree 
of multiprogramming is 3500. For Table IV the total degree of multi- 
programming is 85, but there are 17 classes. As shown, the percentage 
error in both cases is well within an acceptable range. 

In conclusion, the computational experiments suggest that the ap- 
proach based on expansions of integral representations is robust, 
computationally fast, and to be recommended for a variety of problems, 
large and small. 


X. GENERALIZATIONS: INTEGRALS IN NETWORKS WITH MANY 
PROCESSORS 


This section shows that for a large class of closed Markovian 
networks, the partition functions possess simple integral representa- 
tions. This result provides the basis for future work on the computa- 
tions of the integrals from expansions rather like those given in this 
paper. The above-mentioned class of networks allows an arbitrary 
number of service centers with flexibility in operating disciplines. It is 
in fact the same class of networks shown by Baskett et al., in Ref. 4, 
that has the product form in the stationary distributions, except that 
we do not allow the service rate in Type 1 centers to depend on the 
number in queue. (To some extent this is only for convenience because 
for some specific and interesting dependencies, we have obtained the 
integral representations.) 

The representation of the partition function that is obtained is as a 
multiple integral, i.e., as an integral in Euclidean space of dimension q 
where g is the number of queuing centers, which are the centers of 
Type 1, 2, and 4 in the notation of Ref. 4. However, in spite of the 
complexity of the partition function, the form of the integrand is 
remarkable for its simplicity. 
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10.1 Background on product form solutions 


As previously, we let the number of classes of constituents be p. We 
henceforth consistently index the classes by the symbol j; when the 
index for summation or multiplication is omitted, it is understood that 
the missing index is 7, where 1 <7 <= p. A total of s service centers are 
allowed. We will find it natural to distinguish the centers of Types 1, 
2, and 4, which have queuing, from the remaining centers of Type 3, 
which do not. Thus, centers 1 through q will be the queuing centers 
while (q + 1) through s will be the Type 3 centers, which have also 
been called think nodes and infinite server nodes. 

Let the 


stationary state probability = z(yi, yo, ---, Ys), (124) 
Yi A (mii, Nei, ee", Npi), L= i= §, 


n; 4 number of class— jobs in center 1. 


The well-known results on Markovian closed queuing networks with 
product form solutions may be given in the following form:"* 


1 s 
M(Y1, -*+,Ys) = G II zily,), (125) 
i=l 
where wi(Yi) = (inj)! Il () ; l<is qa 
jis 


-1 (2). (q+l<iss. (126) 


nyt 


In the above formulas we have taken into account the previously 
stated assumption, namely, for the first-come-first-served discipline in 
Type 1 centers the service rate is independent of the number of jobs 
in queue. Also, in (126), 


__ expected number of visits of class7 jobs to center i 


yi ’ 


service rate of class 7 jobs in center Z 
where the numerator is obtained from the given routing matrix by 
solving for the eigenvector corresponding to the eigenvalue at 1. 

In (125) G is, of course, the partition function and it is explicitly 


G= Do-s Y TL alya, (127) 


1n,=K, In,=K, i=1 


where we have written 1’n; for }ii-1 nj; and the condition 1’n; = K; to 
indicate the conservation of jobs in each class. Thus, 
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G=))---D fi (« ny) TI FY |, (1 YI. (128) 


10.2 Integral representation 


Using Euler’s integral, see (10), 


20 oo q q Lyi 
-[-[e(-2u) 3 a f(T") | 
0 0 t=] Vn,=K, 1'n,=K, i=1 Nyi- 
| (112) | dn ++ a (129) 


Now by the multinomial theorem, 


vane [fb 


Kj 
-T] (3 Pyilii + ) ei] du, e402. dug. (130) 
i=qt+l1 
It is noteworthy, but not surprising, that the parameters p,; for all the 
Type 3 centers appear lumped together. 

We now introduce the large parameter WN and define 


K; 
B= Wy 1<jsp, (131) 
exactly as in (12). However, we define 
w=(—" |N, 1sjsp, 1sisq (182) 
> Pjm 
m=qtl 


which is reciprocal to the natural extension of the parameters {y;} 
defined in (13). In common with Section 3.1, the suggestion in the 
notation is that in the generic large network {y,;} and {8;} are O(1). A 
tacit assumption being made is that all job classes are routed through 
at least one Type 3 (infinite server) center. 

On substituting (131) and (132) in (180) and after a change of 
variables we obtain a form for the partition function, which is sum- 
marized below. 


Proposition 18: 


e-(fil(SoYm) [Lome om 


634 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1981 


where 


Z= (21, Za, 9°, 2q)’, 
Pp 
f(z) = Vz — > B; log(1 + Tz), 
J=1 


1=[1, 1, coe LY, 
Ty = [ya, Yi2, *%% 5 val’, 1 <j =p. Oj 


As before, the term in brackets in (133) will typically not be required 
to be computed. 

Future work will consider the expansions appropriate for the com- 
putation of the integral in (133). 


APPENDIX A 


Notation for Asymptotic Expansion, Series Reversion 


Asymptotic expansion: A series \%-o a;/N?’ is said to be an asymp- 
totic expansion” of a function F(N) if 


n-1 


F(N) - 2, a;/N/ = 0(N~") as N->o 
for every n = 1, 2, --- . We write 
F(N) ~ ¥ aj/N’. 
j=0 
The series itself may be either convergent or divergent. 


Series reversion: If u = f(z), Uo = f (Zo), f’ (zo) ¥ 0, then by Lagrange’s 
expansion” 


7 ” (u — Uo)’ av Z— 2 : 
panegtsf (ee) 


0 


In particular, if f(-) is specified in a Taylor series, the above expansion 
identifies the coefficients { g;} in the reversed power series z — 2) = 
7-1 B(U — Uo)’. We specifically identify the leading coefficients of two 
reversed series that have been used in Section 4.1 and Section 7.1. 

If 


u=fizt fez? +---, (135) 
then z=gut go +---, (136) 
where 


fig = 1, 
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182 = —fr, 

figs = 2f2 — fifs, 

figs = 5fifeh — fifs — 5f2, 

figs = 6fifefs + 3fifs + 14f2 — fifs — 21f fh. 
Similarly, if 


u = fez? + faz? + faz? + +> (137) 
then z= gil? + gou + guu’? + ++, (138) 
where 
2g = 
fige = “fe 


a5 = (5f3 — 4fofs)/8, 
figs = (8hfsfi — ffs — 2f3)/2, 
43/29, = (—64f3 fe + 224f3 fa fs + 112f3 f2 — 504f2 f3 fi + 2318) /128. 


APPENDIX B 
Proof of Proposition 6 


Before giving the proof, it is worthwhile to generate expressions for 
the leading derivatives of g(-). In the notation of Section 4.1, u = f(z) 
and z= g(u), 


g(u) = 


os 

f(2zV 

E(u) = — f(z)/( F(z)’, 

E(u) = [—fF(z) FP) + 8 Fz)? (FP@))*, 

E(u) = [=F (z)( F(Z)? + 10F(z) F(Z) FZ) 

— 15( fF (z))°1/ (fF (z))’. (139) 


For notational convenience, let y; and ¢; denote, in this appendix only, 
g”(u) and f(z) respectively. Recall that ¢ > 0 and that (—1)‘¢; > 0 
for i = 2, 3, ---. We will show by induction that (—1)"yn < 0, n = 
1, 2,3, +++. 

Let the induction hypothesis be the following 


n even: $2" "yn = +++ — (1)™(X) + (b)27"Y) «++, (140) 


where 0 = 7; max [21, 21 + 1] < 2n — 1; X and Y are arbitrary products 
of $2, $3, +++, $n With arbitrary positive numerical weight and the 
property that X > 0, Y < 0 for all z= 0. 
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n odd: 61" "yn = +++ + (1)"*(U) — (o1)**(V) ++, (141) 


where U and V are like X and Y including that U > 0, V < 0 for all 
z=0. 

For the proof, take the case of n even, first. A key observation 
is that since @, does not occur in either X or Y, dX/dz < 0 and 
dY/dz > 0. That is, in common with the functions ¢2, ¢3, --- the 
derivative has opposite sign from the function. Also, from (140) 


xX Y 


Yn = 8 oe + on are 


(142) 


Differentiate with respect to u, 


1 (= — 2i — 1)d2 9 ) 


aaa gn geri 
1 (-Y(2n—2i-2)d2 YY’ 
+ ry (See + aos epee 
Hence, 
OY" yma = +++ + (h1)"(X (2n — 2i — 1)42) 


— (pi)**1(X’ + Y(2n — 2i — 2)ho) + (i)7**72(Y’) --- (143) 


which has the form appearing in (141) as part of the hypothesis. 
Now take the case of n odd where, from (141), 


Sey WO RN ee (144) 


go ii 
Differentiate and rearrange to obtain, 


Oya = +++ — (b1)"(U(2n — 2i — 1)¢2) 
+ (1)*(U! + V(2n — Bj — 2)h2) — (pi)"7(V’) «es. 


As U’ <0, V’ > 0, the form in (140) of the induction hypothesis is 
satisfied. This concludes the proof. 


APPENDIX C 
Proof of Proposition 10 


In the following we shall need the sign property of the derivatives of 
f(-) as given in (25) and (26). Also recall f, = f’(0) = —a 
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i ze N@dz< | ze Nh? = 1(Na)? (145) 


0 


i) 


fee) 


> ze Nhzth2"\q5 


oo; 


a ae a N/4fe —_ _— 
af \ Nh e (1 — erf(— a VN/4f2)) 
1 
+ INA (146) 
Similarly, 
i e ‘2qz < -1/(aN) (147) 
0 
T e®N/4fo(4 as 
> Vann’ (1 — erf(— a VN/4f2)). (148) 
Thus, 
ze N@dz : 
ditt ae 





[ ody a°N*” 11 — erf(- « VN/4ft)] 
0 


_ 1+ Vit 8f/e?N 
> 2\|a|N : 


where, to bound the term dependent on the error function, we have 
used the left inequality in the following” 


(149) 


(x + Vx? +2) 2 <e* | e dy <[x(1 + V1 + 4/(ax?))]. (150) 


We obtain from (146), (147) and (150) the results analogous to (149): 


i ze Nadz 
0 ee 8 (: 2 ) (151) 


i * oes ofp 1+ V1 + 16f:/(707N) 
0 


Equations (149) and (151) taken together with the representation for 
the utilization given in (22) is the content of Proposition 10. 
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APPENDIX D 
Proof of Proposition 11 


To prove the proposition, we begin by substituting (54) in the 
expression for I in (56) to obtain after a change of variable 


I= (VN +c) | e Nt TT (1 + v/(ai VN))&%t4™ dp, (152) 
0 


We may write (152) as 


I= | e4”/2-Byp(y JN)dp, (153) 
0 


where 
A = ¥bi/a?, B=c — Ydi/ai (154) 
and 
h(v, VN) = (VN + c) exp[Av?/2 + Bp 
— (VN + c) IIT] (1 + v/(ai VN))o%*4],— (155) 


The quantity A(», VN) has been defined so that when the second 
bracketed expression [-] is written as exp(log [-]) and then the log 
[.] term expanded, a cancellation of the leading terms is effected by 
multiplication with the exponential term in (155). In this step, notice 
has to be taken of the constraint in (55) in that exp[Av?/2 + By — 

»(VN + c)] = exp[Av?/2 — v(VN N bi/a; + ¥.di/ai)]. In this manner we 
obtain, 


h(v, VN) = (VN + cyexr( 5 FAS), (156) 
j>1 


where 








1 d; , 
ae 1)/F(v) = (FH a] is ae = Fes +2 2 ai) me, (157) 


In (156) we have an exponential of a power series in 1/VN which may 
equivalently be represented directly as a power series in 1/ VN: 


hiv, VN) = (VN +c) ¥ Gv) /VN%. (158) 
j=0 
For example, Go = 1, Gi = Fi, Go = F?/2 + Fo, G3 = F?/6 + Fi Fo + Fs. 
A feature to observe is that G;(v) is a polynomial of degree 37 in v. 


For the final form, 


ho, WN) = Y Hy /WN", (159) 
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where H,(v) = G;(v) + cG;-1(0), also a polynomial of degree 37 in vp. It 
is understood that G_i(v) = 


Insertion of (159) in (153) yields Proposition 11, namely 


I= y c/ VN, (160) 


J=0 


where a= | e74"/2-BY TT (py) dy, (161) 


0 
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Time-Frequency Multiplexing (TFM) of Two 
NTSC Color TV Signals—Simulation Results 


By B. G. HASKELL 
(Manuscript received September 18, 1980) 


Simple frequency-division multiplexing (FDM) of two television sig- 
nals onto a microwave radio channel is frequently unsatisfactory 
because of nonlinearity-induced crosstalk from one picture into the 
other. With time-frequency multiplexing (TFM), two successive scan 
lines (or fields) of one picture are frequency multiplexed so that they 
can be sent in one line (or field) period. During the next time interval, 
two successive lines (or fields) from the other picture are transmitted, 
thus avoiding crosstalk between pictures. To reduce the bandwidth 
required, one of the two simultaneously transmitted lines (or fields) 
is sent as an analog differential signal. In this paper, we discuss 
computer simulations of these techniques with application to a 20- 
MHz-bandwidth microwave radio channel. Effects of filtering and 
nonlinearities are included insofar as possible. 


1. INTRODUCTION 


It has been known for some time that frequency-division multiplex- 
ing (FDM) of two 4.2-MHz National Television System Committee 
(NTSC) color television signals onto one microwave radio channel is 
often unsatisfactory because nonlinear distortion causes visible 
crosstalk between the two pictures. Also, we suspected that 2:1 time 
compression followed by time multiplexing would lead to problems 
because of the color carrier (initially 3.58 MHz) being transmitted at 
7.16 MHz, where it is much more susceptible to degradations such as 
selective fading. However, there are other techniques which show 
promise of not having such problems. 

Here we propose time-frequency multiplexing (TFM), which involves 
sending two successive scan lines (or fields) from one picture during 
one line (or field) interval, followed by two lines (fields) from the other 
picture in the remaining time interval. The two successive lines (fields) 
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from one picture would be transmitted simultaneously via frequency- 
division multiplexing. Multiplexing in this manner avoids crosstalk 
from one picture into the other which would otherwise be caused by 
nonlinearities in the transmission if simple FpM were used. We hy- 
pothesized that nonlinearity induced crosstalk between successive lines 
(fields) in the same picture will be much less serious because of their 
high similarity and the fact that the ghost of a picture overlayed upon 
itself is practically invisible. With these methods the two incoming 
pictures would have to be in phase synchronism. 

While these methods seem feasible on a qualitative basis, quantita- 
tively many questions arise about their behavior. For example, if one 
of the successive lines (fields) is sent as a differential signal, e.g., line- 
to-line difference or field-to-field difference, how much can the band- 
width of the differential signal be reduced before detectable picture 
degradation results? Also, how much intermodulation crosstalk can be 
tolerated? 

To take a first step toward answering these questions and to gain a 
rudimentary idea as to the limitations of these techniques, we under- 
took computer simulations. The computer facility display had a reso- 
lution of only 256 X 256 pels (picture elements), which is only about 
one-quarter the resolution of an NTSC picture. Thus, each processed 
picture can be regarded as one quadrant of a full-resolution NTSC 
picture. Modulation and filtering inherent in Ntsc color multiplexing 
and demultiplexing were realized by digital filters operating at a 
simulated sampling rate of four times the color subcarrier frequency. 
We maintained full Ntsc bandwidth for luminance and chrominance, 
insofar as possible, even though very few home or studio monitors in 
the U.S. are capable of displaying full resolution (see Appendix A for 
details). We simulated time and frequency multiplexing inherent in 
the transmission methods under discussion using baseband digital 
filters. Effects of nonlinearities were also simulated at baseband, as 
will be described below. 


1. DESCRIPTION OF THE METHOD 


Long-haul broadcast color Tv transmission is presently by FM, as 
shown in Fig. 1. One might suggest transmitting a second vestigial- 
sideband (vsB) signal via FDM as shown in Fig. 2. However, difficulties 
immediately arise: (t) The microwave amplifier nonlinearities can 
cause intermodulation crosstalk, causing ghosts at an unacceptable 
level. (ii) Modulating the baseband plus vsB signal onto an FM carrier 
may exceed the bandwidth allowed, since the modulation is not strictly 
narrowband FM. Carson’s rule estimates the required rf bandwidth in 
Fig. 2 as = 2 (10 MHz + deviation). 

Ghosts caused by nonlinearity-induced intermodulation can be made 
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7fo = 3.58 MEGAHERTZ 
f 


~10 0 4.2 10 
MEGAHERTZ 


Fig. 1—tTVv transmission via narrowband FM. The bandwidth actually used is 8.4 MHz 
plus the bandwidth used by sidelobes (not shown here). If the peak Frm deviation is AF, 
then Carson’s rule estimates the required bandwidth as 8.4 MHz + 2AF. The color 
subcarrier frequency is f. ~ 3.58 MHz. 


much less visible if time-division multiplexing of the two picture signals 
is employed instead of frequency-division multiplexing. As explained 
in the introduction, this entails sending two lines (or fields) from one 
picture in one line (field) interval followed by the same information 
from the other picture in the next time interval. This requires that the 
two incoming pictures be synchronized with each other at all times (if 
this is not the case, a synchronizing device is needed).* Time multi- 
plexing eliminates crosstalk between the two pictures. 

Unfortunately, if FpDM is used to transmit simultaneously two suc- 
cessive lines (or fields), nonlinearity-induced crosstalk between these 
two components will still exist. However, two successive lines (or 
fields) from the same picture tend to be very much alike. Thus, 
overlaying the ghost of one onto the other should be practically 
invisible. 

From Fig. 2, we see that there may not be enough bandwidth to 
send two successive lines (or fields) using simple vestigial-sideband in 
the upper band. However, if a difference signal, e.g., line-to-line differ- 
ence or field-to-field difference, is transmitted in the upper band, then 
the aforementioned difficulties can be reduced. A line (or field) differ- 
ence signal has much less power than the input video and, judging 
from picture coding experience, probably requires less bandwidth for 
transmission as well. It has little or no power near Dc and, therefore, 
lends itself to single-sideband modulation (sss). The arrangement of 
Fig. 3 may be suitable. In this case the differential signal in ssB 
modulated, added to the baseband signal and passed to the radio 
system for transmission via FM. The single-sideband carrier could be 
transmitted during horizontal blanking so that it would not use up 
valuable bandwidth during the visible part of the picture, or it could 
be generated from the color subcarrier. 

In deriving the differential signals, differences should be taken 


* Synchronization of the two pictures should be such that horizontal sync pulses 
align. In this way, switching from one picture to the other can occur, for example, 
following the color burst, and switching transients will not be visible in either picture. 
An ordinary frame synchronizer should suffice. 
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Fig. 2—Transmitting two Tv signals via simple frequency- division multiplexing using 
vestigial-sideband in the upper frequency band. This scheme is unsuitable because of 
(t) nonlinearities and (ii) the fact that the modulation is not strictly narrowband FM and 
may require more than the allowable bandwidth according to Carson’s rule, i.e., FM 
bandwidth ~ 2 x (message bandwidth + deviation). 


between picture elements which on the average have the same color 
subcarrier phase. This prevents a strong color subcarrier component 
from appearing in the difference signal. For example, in Fig. 4 three 
successive lines in one field are shown, and the center one is to be sent 
via a line difference signal. Pels A, B, C, D, and X all have the same 
color subcarrier phase (on a flat-colored area). Thus, a suitable line 
differential signal D, would be X minus an average of A, B, C, and D. 
At the receiver, pels A, B, C, and D are averaged in exactly the same 
way, and the received differential signal is added to recover a replica 
of the original center line. 

In Fig. 5 the center line is from the present field while the other two 
are from the previous field. Pels A, B, C, and X all have the same color 
subcarrier phase on a flat-colored area. A suitable differential signal 
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Fig. 3—Time-division multiplexing (baseband and differential signals are from the 
same source). If one of the two lines (or fields) is sent as a differential signal, then it can 
be placed into the upper band via single-sideband modulation prior to transmission over 
the radio facility. The single-sideband carrier could be transmitted sometime during 
horizontal blanking so that it would not use up valuable bandwidth during the visible 
part of the picture, or it could be generated from the color subcarrier. 
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Fig. 4—Three successive lines in one field are shown, and the center one is to be sent 
via a. line-difference signal. Pels A, B, C, D, and X all have the same color subcarrier 
phase (on a flat colored area). Thus, a suitable line differential signal Di would be X 
minus an average of A, B, C, and D,ie., Dp = X — 4(A+B+C+D). 


Dr is X minus an average of A, B, and C. Since C is closer to X than 
A or B, it should be weighted relatively more, i.e., a > % should be 
used. At the receiver, we recover the center line by forming the 
weighted average and adding the received differential signal. In both 
Figs. 4 and 5 the pel spacing implies a sampling rate of four times the 
color subcarrier frequency, which is the rate used in the computer 
simulations. 


ili. BANDWIDTH REQUIRED FOR DIFFERENTIAL SIGNALS 


The first question which we confronted in the computer simulations 
was “How much can the differential signal be bandlimited before 


Fig. 5—Three successive lines in one frame are shown, and the center one is to be 
sent via a field difference signal. The center line is from the present field while the other 
two are from the previous field. Pels A, B, and C all have the same color subcarrier 
phase on a flat area. A suitable differential signal Dry is X minus an average of A, B, and 
C, ie., Dr = X — [aC + (A + B)(1 — a)/2]. Since C is closer to X than A or B, it should 
be weighted relatively more, i.e., a > % should be used. At the receiver the center line 
is recovered by forming the weighted average and adding the received differential signal. 
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detectable picture distortion results?” In the simulations, we main- 
tained full NTsc resolution in the color carrier modulation and de- 
modulation (see Appendix A for details), and used a monitor with 
separate red, blue, green (RBG) inputs for display. We used slides 
from the Society of Motion Picture Television Engineers (SMPTE) with 
R, B, and G components digitized to eight bits by means of a flying 
spot scanner for the pictures. 

Starting with RBG signals for each picture, we first produced NTSC 
composite signals as described in Appendix A. From these, we obtained 
reference pictures by NTSC demodulation into the R, B, G components 
required for display. 

To test line-differential transmission, we computed the line differ- 
ence shown in Fig. 4, from the NTsc composite signals for alternate 
lines in each field. We then bandlimited the differential signal to 
various bandwidths, using finite impulse-response digital filters to 
avoid phase distortions. The bandlimited differential signals were then 
added to the appropriate averaged signals in alternate lines to recover 
replicas of the original baseband composite signals. The composite 
signals were then demodulated to obtain processed pictures. 

We compared each processed picture with its corresponding refer- 
ence picture by switching between the two on the same monitor and 
closely inspecting small areas for perceivable change. From this we 
concluded that as long as the line differential signal had a bandwidth 
of at least 3 MHz, no perceivable difference would result. Below 3 
MHz, a slight vertical color shading could be detected in a few areas, 
such as the boundary beween red lipstick and white teeth. | 

To test field-differential transmission, we carried out the same 
procedure using the differential signal shown in Fig. 5. We examined 
several values of a in the range 0.4 to 0.6, but the results were the 
same in all cases. We concluded that as long as the field differential 
signal had a bandwidth of at least 2 MHz, no perceivable difference 
would result between the reference pictures and the processed pictures. 
To test the effect on movement rendition, we carried out the processing 
on sequences of frames moving in pure translation, vertically and 
horizontally. The results were the same as for still frames and are 
summarized in Table I. 

It is certainly possible to generate test signals which cannot be 
transmitted without degradation by the methods suggested here. For 
example, monochrome vertical stripes at the color subcarrier frequency 
generate a strong line-differential signal at the color subcarrier fre- 
quency which would not pass through a 3-MHz low-pass filter. If the 
vertical interval test signal (vITs) happens to fall on a line that is 
transmitted via a bandlimited differential signal, then degradation 
results. The vitTs is, after all, designed to test baseband transmission 


648 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1981 


Table I—Bandwidths required 
for transmitting differential 
signals if no perceivable picture 
degradation is allowed. 


: Required 
Signal Bandwidth 
Baseband Video 4 MHz 
Line Difference 3 MHz 
Field Difference 2 MHz 


and horizontal resolution. As such, it should be sent intact without 
any differential processing. Other test signals should be devised to 
measure vertical color resolution in systems, such as proposed here, 
which limit the transmitted vertical color resolution to that actually 
required by pictures. 


iV. APPLICATION TO A 20-MHz BANDWIDTH MICROWAVE RADIO 
CHANNEL 


4.1 Channel characteristics 


Figure 6 shows a typical long-haul, microwave, television transmis- 
sion system. Preemphasis of high frequencies gives some protection 
against nonlinearities. Figure 7 gives the preemphasis characteristic; 
de-emphasis is the inverse of preemphasis. The FMT produces 32.0- 
MHz deviation per volt input. The gain a is chosen so that a nominal 
deviation of about 2 MHz results from a one-volt peak-to-peak color 
bar video signal. However, sustained deviations of 3.2 MHz (found by 
computer search) can occur for certain input signals, and transient 
deviations as high as 3.7 MHz (preemphasized 2T pulse)’ are possible. 
Ac coupling occurs before the FMT. However, this has little effect on 
possible peak deviations since its passband extends well below most 
frequencies of interest. 


H(0)=b H(co)=1 


FMT 
VIDEO 
IN FM 
PREEMPHASIS MODULATOR 
TRANSMISSION 


FMR 
VIDEO 
OUT GAIN y FM 
d DEMODULATOR 


Fig. 6—A typical long-haul television transmission system. a = —12.6 dB, b = —13.3 
dB, c = 16 dB, d = —3.4 dB, acd = 1. De-emphasis is the inverse of preemphasis. The 
FM modulator produces 32.0-MHz deviation per volt input. 
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Fig. 7—Preemphasis characteristic. H(s) = (w: + s)/(we + s), where f; = 199 kHz, 
fz = 873 kHz, and H(0) = 6 = —13.3 dB. Preemphasis of high frequencies gives some 
protection against nonlinearities and improves SNR in colored areas of the picture. 


Nonlinearities vary considerably between channels. However, one 
such microwave radio system (4000-mile terrestrial, 4GHz) was mea- 
sured, and Appendix B estimates video degradation for a standard 
vits. Differential gain is estimated at 8.5 IRE* (13 IRE allowed), and 
worst case differential phase is estimated at 9.7° (5° allowed). These 
results indicate that the measured channel was marginal and, there- 
fore, simulation results based on these measurements should be on the 
conservative side. 


4.2 Time-frequency multiplexing 


We wish to transmit a baseband signal of 4.2-MHz bandwidth and 
a single-sideband modulated differential signal in the positive fre- 
quency band as shown in Fig. 3. Figure 8 shows a configuration for 
doing so. It is similar to Fig. 6 except that a 20-MHz bandpass filter 
following the FM modulator is shown explicitly, and the preemphasis 
is changed to admit the possibility of attenuating the higher frequen- 
cies of the baseband signal somewhat to obtain satisfactory perform- 
ance. 

A line differential signal, which requires 3-MHz bandwidth, is prob- 
ably not suitable for this application. However, a field differential 
signal, which requires only 2 MHz of bandwidth, could fit into the 5 to 
7-MHz band if sharp cutoff filters were used. Hereafter, results will 
apply to a system using field differential signals. 

The differential signal has no dc component and is small most of the 
time with ordinary pictures, occasionally reaching values of +0.36M.t 


* 1 volt = 140 IRE. 
+ Measured by computer simulation using SMPTE slides. Electronically generated 


tests signals can be concocted, however, which produce differential signal values as large 
as +M. Here, M = 0.714 volt (see Fig. 9). 
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Carson’s rule states that the sum of the total message bandwidth and 
the rm frequency deviation should not exceed one-half the available 
channel bandwidth. Since the channel bandwidth is 20 MHz, and the 
proposed message bandwidth is 7 MHz, the rm deviation should not 
exceed 3 MHz. Unfortunately, as we have seen, the system in its 
present form already has deviations larger than 3 MHz, even without 
the addition of a modulated differential signal. However, these devia- 
tions are often short lived and occur at sharp edges in the picture 
where considerable distortion can be withstood. Still, some signal 
attenuation is obviously necessary to accommodate a 7-MHz message 
bandwidth. However, it may not be necessary to attenuate the signals 
so.much that the deviation is always below 3 MHz. If the FmrT is 
followed by a reasonably good 20-MHz bandpass filter, then some of 
the responsibility for maintaining bandwidth might be shifted to it, at 
the cost of some picture degradation. This degradation due to the 
bandpass filter will only occur if the preemphasized baseband signal is 
large, while at the same time the differential signal is large. In this 
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Fig. 8—Configuration for transmitting a single-sideband modulated differential signal. 
a = —12.6 dB, 6 = —13.3 dB, c = 16 dB, d = —3.4 dB, acd = 1. It is similar to Fig. 6 
except that a 20-MHz bandpass filter following the Fm modulator is shown explicitly, 
and the preemphasis is changed to admit the possibility of attenuating the higher 
frequencies of the baseband signal somewhat to obtain satisfactory performance. Gain 
e causes attenuation of the modulated differential signal (including the effect of single- 
sideband modulation), and gain f causes attenuation of the high frequencies of the 
baseband signal. e = f = —3 dB gives just visible distortion. 
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case, and in this case only, the bandpass filter will attenuate the 
modulated differential signal, thus causing distortion. This distortion 
should be relatively invisible, however, due to its proximity to sharp 
brightness transitions in the picture. In Fig. 8, gain e causes attenuation 
of the modulated differential signal (including the effect of single- 
sideband modulation), gain f causes attenuation of the high frequencies 
of the baseband signal. Suitable values for e and f must be determined 
experimentally. However, some guidelines can be obtained intuitively 
as follows. 

If a constant peak-white baseband value of M = 0.714 volt (=100 
IRE) occurs at the same time as a differential signal peak of 0.36M, 
then the input to the FM modulator (assuming worst case Dc = 0)* is 


x = abM + 0.36aM cos’wpt, (1) 


where wp is the differential signal carrier frequency, and cos’ indicates 
single-sideband modulation. If the conservative position is taken, of 
never allowing a deviation above 3 MHz, then from values previously 
defined, e should be less than —0.4 dB. This is a very modest atten- 
uation. 

Let us consider another input signal. A maximum sustained preem- 
phasized baseband signal occurs if the video input is 0.52M + 0.48M 
cosw,t. If this value occurs at the same time as a peak differential 
signal value of 0.36M, then again assuming the worst case (Dc = 0), 


x = 0.52abM + 0.48afM cosw-t + 0.36aeM cos’wpt. (2) 


If, in this case, the very conservative approach of 3-MHz maximum 
deviation is taken, then with an FMT deviation of 32 MHz/volt, 


0.52abM + 0.48afM + 0.36aeM = 0.094 volt. (3) 
If equal values are chosen for e and f, 
e= f= —5.4 dB. (4) 


Although choosing e and f to satisfy Eq. (4) eliminates almost all 
distortion caused by the 20-MHz bandpass filter, it also leads to a 
decrease in overall SNR, i.e., fade margin. Since some of the bandpass 
filter distortion will be invisible, it seems wasteful to use such low 
values for e and f. Larger values should be feasible with the same 
subjective picture quality. Thus, another of the reasons for performing 
the computer simulations was to estimate the range of acceptable 
values for e and f. 


* Average signal level. This is not transmitted in an ac coupled system. 
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V. COMPUTER SIMULATION OF NONLINEARITIES 


Since the transmission channel is not precisely understood and is 
time varying as well, such simulations can only be expected to give a 
vague idea of overall performance. Hopefully, however an estimate of 
worst-case performance can be obtained if a conservative approach is 
taken. The channel characteristics used in Appendix B and shown in 
Table II probably represent a marginally acceptable situation insofar 
as nonlinearities are concerned. Thus, the use of these data in the 
simulations should give conservative results. 

There are several approaches to modeling the nonlinearity, including 
the rather complicated procedure of numerically carrying out the 
complete ssB modulation, multiplexing and FM operations, and passing 
the result through the nonlinearity. However, the simulations to be 
described below were all carried out at baseband for simplicity and 
since intermodulation crosstalk was the main item of interest. 

The basic approach taken in evaluating the effects of nonlinearities 
was to introduce degradations into the NTSc composite signal and see 
if they were visible in the resulting reproduced picture. Two types of 
degradation were studied: (z) intermodulation crosstalk resulting from 
nonlinearities, and (it) the effect of the 20-MHz bandpass filter follow- 
ing the FM modulator in Fig. 8. 

Baseband and differential signals were first computed as in Section 
III for each of the original SMPTE pictures. Intermodulation crosstalk 
terms were then estimated as described in Appendix C. Some of the 
terms degrade the baseband signal, and some degrade the differential 
signal. 

The effect of the 20-MHz bandpass filter was approximated as 
follows: If at certain times the sum of the baseband and differential 
signal magnitudes at the input to the FMT was large enough to produce 
a 3-MHz deviation (again assuming worst case dc = 0), then the 
magnitude of the differential signal was reduced until the deviation 
fell below 3 MHz. If the baseband signal by itself was large enough to 


Table Il—Results of two-tone measurements on a typical long- 
haul network. A = B = —26 dBv, c = 16 dB. From the 
measured harmonic outputs, the coefficients u and v of the 
third-order model can be obtained. 


ucAB FMR Vik ewan 


average average dBv 
a(MHz) (MHz) dBv at 8 u v 
Bt+aand at B + 20: 
poo and 8 — 2a 
0.1 6.55 —41 —49 0.56 5.96 
0.1 3.55 —35 —46.5 1.12 7.94 
0.1 1.55 —36 —39 1.00 18.84 
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produce a 3-MHz deviation, then the differential signal was set to 
zero.” 

In the simulations, gain factors a through f were all incorporated as 
shown in Fig. 8; factors a through d were fixed having the values 
shown, while e and f were made variable. The nonlinearity coefficients 
were taken from Table II, depending on the applicable frequency 
range. 

Results were generally encouraging, considering the pessimistic and 
conservative assumptions built into the simulations. Using attentua- 
tion factors e = f = 0.5 resulted in no visible degradation, as expected 
from eq. (4). Larger values in the range 0.7 to 0.8 (—3 to —2 dB) gave 
rise to barely visible color shading in a few areas such as the boundary 
between red lipstick and white teeth. Otherwise, the effects were 
invisible. 


Vi. CONCLUSION 


Color pictures were processed via computer simulation to test the 
feasiblity of using time-frequency multiplexing to transmit two broad- 
cast-quality color television signals over a 20-MHz bandwidth micro- 
wave radio channel. We considered bandwidth allocation and nonlin- 
earities in detail. At every turn in the study, conservative assumptions 
were made. The displayed pictures produced full-color bandwidth even 
though very few home or studio monitors in the U.S. are capable of 
doing so. Nonlinearity parameters which were used in the simulations 
were obtained from an, at best, marginal channel. Worst-case theoret- 
ical nonlinearity impairments were assumed. And, finally, for picture 
quality assessment an A-B comparison on the same monitor was used, 
which is one of the most stringent tests known. 

Results were rather encouraging. There appeared to be enough 
bandwidth to accomplish the transmission. With a modest 3-dB re- 
duction of signal level, simulated transmission defects were practically 
invisible. 

Utilization of field differential signals requires field memories at the 
transmitter and receiver. Although they are expensive at present, their 
cost is dropping. Moreover, if synchronization of two pictures is re- 
quired at the transmitter, then the memory of the synchronizer can be 
used. 

Audio transmission has not been considered here. One might be 
tempted to place audio carriers in the vacant frequency band between 
the baseband and differential signals. However, this would require a 
rather linear channel, which is the very problem that the techniques 


* An automatic gain control could, in fact, be implemented in this way to guarantee 
deviation less than 3 MHz. 
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of this paper are meant to avoid. A better solution would be to transmit 
the audio digitally, either in the horizontal sync pulse? or during 
vertical blanking. 


APPENDIX A 
Computer Simulation of NTSC Color Multiplexing and Demultiplexing 


The computer simulation facility has a picture resolution of 256 x 
256 pels per frame. At a frame rate of 30 Hz, 15 lines vertical blanking, 
17 percent horizontal blanking, this implies a sampling rate of ~2.51 
MHz. Thus, video signal bandwidth of the computer is conservatively 
estimated to be about 1 MHz. Since U. S. broadcast rate NTsc color 
television signals have a nominal bandwidth 4.2 MHz, scaling by a 
factor of approximately four is required to do meaningful simulations. 

The ntsc broadcast color subcarrier frequency is approximately 3.58 
MHz. Thus, for the simulations, a value around 0.89 MHz should be 
used. In addition, the simulation color subcarrier frequency must be 
an odd multiple of half the line scanning frequency of the simulation 
system (fx ~ 8.1 kHz). Thus, a convenient value for color carrier 
frequency that was chosen for the simulations was f, = 108.5fy ~ 0.88 
MHz. 

The simulated sampling rate of 4 f. equals 434 fy or 434 pels per line. 
Assuming about 17 percent horizontal blanking, we have 358 pels in 
the active area of the picture (which corresponds to the 256 pels of the 
display). To convert the original red, blue, and green (RBG) signals 
having 256 pels per line to the desired 358 pels per line, linear 
interpolation was used. This sampling rate conversion is suboptimal in 
that a frequency rolloff of 4.7 dB at 1 MHz is introduced. The extra 
effort required to design optimal digital filters for sampling rate con- 
version did not seem worthwhile at this time. After processing, the 
conversion back to 256 pels per line was again performed via linear 
interpolation. 

From the RBG signals, the full bandwidth Y, J, Q signals* were 
obtained by the standard matrix relation® 


Y 0.3 O11 059] R 
I}=|06 -032 -0.28|] B (5) 
Q 0.21 031 -052] G 


The I and Q signals were then bandlimited to ~0.375 MHz and 0.125 
MHz, respectively (0.25 X NTSC values) using symmetrical finite-im- 
pulse-response digital filters to avoid phase distortion. Filters were 
chosen to meet NTSC specifications, yet to be short enough to avoid 
excessive ringing. 


* Y = luminance, J and Q = chrominance signals. 
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Samples at 4 f. were assumed to occur at the peaks and zeros of the 
I and Q carrier sinusoids. The chrominance (CH) pels could then be 
formed simply from 


+CH(N) = I(N) cos(Na/2) 
+ Q(N) sin(N7/2), N=1,---, 358, (6) 


the plus and minus signs being taken on alternate lines of the field. 
The chrominance and luminance were then added together and filtered 
to a bandwidth of 1 MHz, again with a FIR digital filter to avoid phase 
distortion. The result is the composite NTsc color signal. 

Recovery of the R, B, G signals from the composite signal is not 
quite so straightforward if full bandwidth for the luminance and 
chrominance is to be maintained. Comb filtering of three successive 
lines (—0.25L, + 0.5L2 — 0.25L3) followed by bandpass filtering (center 
frequency f.) was used to recover the chrominance signal for the center 
line. Simple subtraction of the chrominance from the composite signal 
yielded the luminance pels. 

Recovery of Jand Q from the chrominance started with a quadrature 
AM demodulation, 


I(N) = +CH(N) cos(Nz/2), 
Q(N) = +CH(N) sin(N2/2), (7) 


the plus and minus signs being taken on alternate lines. The required 
low-pass filtering was implemented simply by replacing zero samples 
by an average of their neighbors. 

At this point, the J signal required equalization since it lost part of 
its upper sideband during the 1-MHz low-pass filtering of the composite 
signal. The frequency range from 0.125 MHz to 0.375 MHz must be 
increased by 6 dB. This was accomplished, again, by a FIR digital filter. 
Having obtained Y, J, and Q, the R, B, G signals were produced by 
inverting the matrix equation (5) above. 


APPENDIX B 
Effects of Nonlinearities on the Vertical Interval Test Signal (VITS) 
In Fig. 6, nonlinearities are often modeled as third-order polyno- 
mials,‘ 
y = f(x) = e(x — ux? — vx’). (8) 


This model is reasonable for transmission systems with “small” non- 
linearities. The nonlinearity coefficients u and v are often estimated 
by sending the two-tone signal 


x =A cosat + B cosft, (9) 
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in which case y/c consists of a number of sum and difference frequen- 
cies as follows:* 


pc: u(A? + B’)/2 (deleted at FMR) (10) 
First Order: A cosat + B cost + very small terms, (11) 
Second Order: —%u(A? cos2at + B* cos2Bt) 
—uAB[cos(a + B£)t + cos(a — B)é], (12) 
Third Order: —%4v(A* cos3at + B® cos3Bt) 
—%vA*B[cos(2a + B)t + cos(2a — B)t] 
—%vB*A[cos(2B + a)t + cos(2B — a)é]. (13) 


In 1974, two-tone measurements were made of a typical long-haul 
network. From these measurements the coefficients u and vu can be 
estimated. Results are frequency dependent and are shown in Table II 
for A = B = —26 dBv,* and c = 16 dB. 

These data can be used to estimate degradations which would occur 
in the transmission of a vertical interval test signal (viTs). Let M = 
100 IRE = 0.714 volt be the blanking-to-peak magnitude of the video 
signal as shown in Fig. 9. Values of gains a, 6, and c are given in Fig. 
6. 


Differential gain (Ref. 1, Section 3.13) 
dc = 0, 
Peak input = M(0.9 + 0.2 cosw-é). (14) 
Thus, 
x = aM(0.9b + 0.2 coswt) 
SA+B coset, (15) 


where A is the low-frequency magnitude, and B is the color-subcarrier 
magnitude. From eqs. (11), (12), and (13) with a = 0 and B = «,, the 
color-subcarrier terms in y/c (when phases align,” which is worst case) 
are 


B — 2uAB — %vA’B. (16) 

Thus, at w, the signal and distortion components in y are, respectively, 
R= Be = 0.2aMc, 

AR = —2uABc — %vA*Be. (17) 


* Zero to peak. 
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PEAK _ 
LUMINANCE ~ 





SYNC-~ 


Fig. 9—Composite video signal, including sync, color burst and, video. S = —0.286 
volt. According to specifications, M can be as large as 0.857 volt. However, to produce, 
the normally used one-volt peak-to-peak composite signal, M = 0.714 volt (=100 IRE). 


De-emphasis has no effect at this frequency. Thus, according to the 

test procedure we should divide the above equations by 0.2ac to bring 

the signal component up to M. Then, peak differential gain is given by 
—2uABc — *vA*Be 


Using values defined previously and u = 1.12, v = 7.94 from Table II, 
AG = —0.06 volts = —8.5 IRE, (19) 
which is well within the 15 IRE allowed.’ 


Differential phase (Ref. 1, Section 3.14) 


The same input signal is used as in the differential gain test. Using 
the same definitions for signal R and distortion AR, worst-case differ- 
ential phase occurs when these two components are exactly 90° out of 
phase.* In this case, peak-to-peak differential phase is 


AR 3 
Ad =2 tan”! R =2 tan”! A(z + 9 vA), (20) 


which for the values defined previously evaluates to 
Ad = 9.7°, 


which is somewhat outside the 5° allowed. However, this is a worst- 
case maximum and might not occur in practice. 


Conclusion 


From these estimates of distortion one concludes that the channel 
measured was not of extremely high quality and may have been only 
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of marginal quality. Thus, simulation results using those measurements 
should be on the conservative side. 


APPENDIX C 
Important intermodulation Crosstalk Terms 


First consider crosstalk between baseband color-subcarrier fre- 
quency and the differential signal. Let BB. be the input baseband 
chrominance component (frequencies = w, ~ 3.58 MHz) and let D be 
the input differential signal (carrier frequency wp ~ 7 MHz). Then, 
from Fig. 8, 


x = af- BB. cosw,t + aeD coswpt 
A A; cosw.t + B coswpt, (21) 


as ‘in eq. (9). Then, from eqs. (11) to (13) the corresponding in-band 
distortion components in y are 


—%ycA? cos2w,t — uAiB cos(wp — w)t 


— %vcA3B cos(wp — 2w,)t. (22) 


Combining eqs. (21) and (22) yields estimates of interfrequency 
crosstalk caused by the above frequency components. In the computer 
simulations BB, was obtained by first bandpass filtering the composite 
signal to retain only components from 2 to 4 MHz. The result was 
then shifted down in frequency by 2 MHz to obtain the signal BB, 
used in eqs. (21) and (22) to compute distortion. 

Similar crosstalk occurs between baseband low frequencies and the 
modulated differential signal. However, because of preemphasis, A in 
eqs. (11) to (18) is small, making this distortion negligible. 

Crosstalk between baseband midfrequencies, e.g., wy around 1 MHz, 
and the modulated differential signal can be estimated similarly. Let 
BBwm be the input baseband component in the midrange frequencies 
(wa = 1 MHz), and define D as above. Then from Fig. 8, 


x = af-BBm coswyt + aeD coswpt 
A Ao coswut + B coswpt. (23) 
From eqs. (11) to (18) the in-band distortion components in y are 
—'%4yucA3 cos2wyt — ucA2B cos(wp — wm)t 
—Y¥,cvA3 cos8wyt — %vcA3B cos(wp — 2wy)t. (24) 


In the computer simulations, BBy was obtained by low-pass filtering 
(2-MHz cutoff) the composite signal. 
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Overflow Models for Dimension® PBX Feature 
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We present numerical results for some traffic overflow systems with 
queuing. Traffic is offered by two independent streams to two groups 
of trunks, with a finite number of waiting spaces for each, and some 
overflow capability from the primary group to the secondary group. 
We consider three different overflow systems, two of which model 
feature packages offered in Dimension® psx. For given offered loads 
and unit mean holding time, we determine the number of trunks and 
waiting spaces in the two groups so that the blocking probabilities, 
and the average delays of queued calls do not exceed prescribed 
values. We also calculate various other quantities, such as the occu- 
pancies of the trunk groups and the probability of overflow from the 
primary group to the secondary group. Finally, we examine the effect 
on the blocking probabilities and the average delays of varying the 
loads offered to a given system. 


1. INTRODUCTION 


In this paper we present numerical results for some traffic overflow 
systems with queuing. Traffic is offered by two independent streams 
to two groups of trunks with a finite number of waiting spaces for each, 
and some overflow capability from the primary group to the secondary 
group. The holding times of the calls are independent, and exponen- 
tially distributed. In two of the three systems considered, the overflow 
capability models feature packages (FP) offered in Dimension® PBx. 
The third system, which is considered for comparison, differs from the 
other two systems in that no overflow is permitted if there is a waiting 
space available in the primary queue. 

Since there are a finite number of trunks and waiting spaces in each 
group, arriving calls may be blocked and cleared from the system. For 
given offered loads and unit mean holding time, we determine the 
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number of trunks and waiting spaces in the two groups so that the 
blocking probabilities and the average delays of queued calls do not 
exceed prescribed values. We also calculate various other quantities, 
such as the occupancies of the trunk groups and the probability of 
overflow from the primary group to the secondary group. In addition, 
we examine the effect on the blocking probabilities and the average 
delays of varying the loads offered to a given system. 

The numerical results are based on different techniques developed 
by Kaufman! and Morrison.”* The basic problem is to solve a large 
sparse system of linear equations for the steady-state probabilities of 
the number of calls in the two groups. Kaufman used a numerical 
technique involving matrix separability (block diagonalization), and in 
addition obtained numerical solutions by means of successive over- 
relaxation techniques. Kaufman has also applied her techniques to 
overflow systems with more than two groups. Morrison confined his 
attention to overflow systems with two groups, and he adopted an 
analytical approach which considerably reduces the dimensions of the 
problem. 

We have been able to obtain very accurate numerical results by the 
methods presented in this paper. Indeed the various steady-state 
quantities of interest obtained by the procedures of Kaufman’ and 
Morrison?” agree to many significant figures. These results are consid- 
erably more accurate, and less expensive to obtain, than simulation 
results. They may be used as a check on the accuracy of approximate 
methods which have been developed for dealing with overflow prob- 
lems, some of which are mentioned later in this section. 

We now describe in detail the overflow systems which we consider, 
as depicted in Fig. 1. There are 7; trunks and q, waiting spaces in the 
primary group, and nz trunks and qg2 waiting spaces in the secondary 
group. Traffic is offered to the two groups by two independent Poisson 
streams S, and S», with rates Ai and Az, respectively. The holding times 
of the calls are independent, and exponentially distributed with mean 
1/p. If all n2 trunks in the secondary group are busy when a call from 
stream S2 arrives, the call is queued if one of the g2 waiting spaces is 
available, otherwise it is blocked and cleared from the system. Calls 
waiting in the secondary queue are placed in service on a first-in first- 
out basis as secondary trunks become available. 

Three cases are considered for the treatment of calls offered to the 
primary group. The first two cases model feature packages offered in 
Dimension PBx. For these two cases, if all m; trunks in the primary 
group are busy when a call in S; arrives, it is placed in service in the 
secondary group if there is a trunk available and there are no calls 
waiting in the secondary queue. If no trunk is available, then the call 
is queued in the primary group if one of the gi waiting spaces is 
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Fig. 1—Mean flow rates for an overflow system with queuing. 
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available, otherwise it is blocked and cleared from the system. The 
treatment of calls waiting in the primary queue depends on which 
feature package is being modeled. Corresponding to FP8, calls that are 
waiting in the primary queue are not allowed to overflow to the 
secondary group, but must wait for a trunk in the primary group to 
become available. Corresponding to Fp4 (or FP7), calls waiting in the 
primary queue may be served by an idle trunk in the primary group, 
or by an idle trunk in the secondary group, provided that there are no 
calls waiting in the secondary queue, which have priority. The overflow 
systems corresponding to FP8 and Fp4 are depicted schematically in 
Figs. 2 and 3, respectively. 

The two cases considered above, although an idealization of the 
actual situation, embody the essential features of the packages. For 
comparison, we consider a third case, which we denote by FPO and is 
depicted schematically in Fig. 4. In this case, if all m, trunks in the 
primary group are busy when a call in S; arrives, the call is placed in 
the primary queue if one of the g; waiting spaces is available. As with 
FP8, a call that is queued in the primary must wait for a trunk in the 
primary group to become available. If all m; trunks in the primary 
group are busy and all q; waiting spaces are occupied when a call in S; 
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Fig. 2—Schematic of the overflow system corresponding to FP8. 


arrives, it is placed in service in the secondary if there is a trunk 
available and there are no calls waiting in the secondary queue, 
otherwise it is blocked and cleared from the system. Note that no 
overflow is permitted if there is a waiting space available in the 
primary queue. This restriction was invoked by Anderson.‘ 

We remark that the systems corresponding to FP4 and FPO are 
particular cases of the system considered by Rath in connection with 
ACD-Ess (automatic call distributor-Ess).° He considered a system 
composed of two queues,.in which one of the queues is allowed to 
overflow to the other, under specified conditions involving the queue 
lengths. He obtained numerical solutions in the case corresponding to 
Fp4 by using a Gauss-Seidel iteration technique. He also developed an 
approximate procedure for analyzing his system based on the use of 
the interrupted Poisson process (IPP).° An approximate analysis of the 
system corresponding to FP4 has been given by Crater under the 
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assumption that the number of waiting spaces in each queue is unlim- 
ited.’ More recently, Shulman described a method of iteration and 
successive approximation, using the IpP as a traffic model, for analyzing 
the system corresponding to FP8 with several bands of groups.® He has 
used this method in the optimal design of facilities for Dimension PBx 
with overflow and queuing features. 

Section II outlines the numerical procedures for evaluating the 
quantities of interest, based on the analytical results of Morrison.”” In 
Section III two iterative techniques used by Kaufman to obtain nu- 
merical results are discussed.’ The numerical results are discussed in 
Section IV. These results were obtained by Morrison’s method and at 
least one of Kaufman’s two methods. 


ll. AN ANALYTICAL PROCEDURE 


We begin by outlining the numerical procedures for evaluating the 
quantities of interest based on the analytical results of Morrison.”” 
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Fig. 3—Schematic of the overflow system corresponding to FP4. 
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Fig. 4—Schematic of the overflow system designated as FPO. 


The interested reader may consult these papers for the details of the 
analysis. Let p;; denote the steady-state probability that there are 7 
calls in the primary and / calls in the secondary, either in service or 
waiting. These probabilities satisfy a set of generalized birth-and-death 
equations, which take the form of partial difference equations con- 
necting nearest-neighboring states. The basic technique is to separate 
variables in regions away from certain boundaries of the state space. 
This leads to some eigenvalue problems for the separation constant. 
The eigenvalues are roots of polynomial equations, and they may be 
evaluated numerically with the help of some interlacing properties.’ . 
The probabilities p;; are then represented in terms of the corresponding 
eigenfunctions, which can be calculated from simple recurrence rela- 
tions. The constant coefficients in these representations are deter- 
mined from the boundary conditions (one of which is redundant), and 
the normalization condition that the sum of the probabilities is unity. 
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When engineering the system, there are various steady-state quan- 
tities of interest that can be expressed in terms of the probabilities p,;. 
The quantities include the blocking (or loss) probabilities L; and La, 
and the average delays (or mean waiting times), W, and W,, of calls 
that enter the queues. These quantities may be expressed directly in 
terms of the constant coefficients that occur in the representations for 
the probabilities p,;;. Thus the steady-state quantities of interest can 
be calculated directly, once the coefficients have been determined from 
the boundary and normalization conditions, without having to calcu- 
late the probabilities pi;. 

For FPO there is just one set of ne + go + 1 eigenvalues and 
corresponding eigenfunctions. There are then nz + g2 + 1 constants to 
be determined numerically from the boundary and normalization 
conditions. We remark that in this case there are (m + qi + 1)(m2 + 
qg2 + 1) probabilities p;;, so that the analytical approach considerably 
reduces the dimensions of the problem. For Fp4 there is an additional 
set of qi eigenvalues, and corresponding eigenfunctions. In this case 
there are gi + n2 + ge + 1 constants to be determined numerically 
from the boundary and normalization conditions. This compares with 
the gi(q2 + 1) + (mi + 1)(m2 + ge + 1) probabilities p,;;. We note that 
there are fewer nonzero probabilities p;; for FP4 since it is impossible 
for calls to be waiting in the primary queue when there is an idle trunk 
in the secondary and no calls are waiting in the secondary queue. For 
FP8 this is not the case, and the probabilities p;; in the corresponding 
region of state space are expressed in terms of gi additional constants 
and a fundamental solution of a partial difference equation. However, 
it is possible to solve for these additional constants in terms of the 
other gi + 22 + gz + 1 constants, by inversion of a triangular matrix. 

Some conservation relations were used as a check on the accuracy 
of the numerical calculations. These involve the steady-state quantities 
depicted in Fig. 1. Thus @; and Q, are the probabilities that calls from 
streams S; and S», respectively, are queued upon arrival. The mean 
departure rate from the primary queue to the primary trunks is Ri, 
and the mean rate of overflow from the primary queue to the secondary 
trunks is Ry. Note that Riz = 0 for both FPO and FP8. The mean 
departure rate from the secondary queue, to the secondary trunks, is 
R22. Since the mean rate of arrival of calls at each queue is equal to 
the mean departure rate, we have 


AQ, = Ru + Ry, A2Q2 = Ro. (1) 


Also, he is the probability that a call from stream S; overflows 
immediately. Moreover, let Xi and X2 denote the average number of 
calls in service in the primary and secondary groups, respectively. 
According to Little’s formula, the average number in a queuing system 
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is equal to the average arrival rate to that system times the average 
time spent in that system.” If we apply this result to the primary and 
secondary trunk groups, it follows that 


pX1 =Ai(1 — Li — Te) — Ry, pX2 = A2x(1 — Le) + Aitiez+ Riz, (2) 


since 1/p is the mean holding time. Since the various quantities in (1) 
and (2) may be expressed in terms of the constant coefficients that 
occur in the representations for the probabilities p;;,”* these relation- 
ships provide a useful numerical check. 

A package of computer programs has been developed by Seery to 
calculate the various steady-state quantities of interest, based on the 
procedures outlined above. The documentation provides the informa- 
tion necessary to use the programs.’° 


lll, SPARSE MATRIX TECHNIQUES 


The birth-and-death equations which determine the steady-state 
probabilities, and which are mentioned in the beginning of Section II, 
may be written as 


Ap = 0, (3) 


where A is a singular, nonsymmetric matrix with negative off-diagonal 
elements and column sums equal to zero. It has only five nonzero 
diagonals. 

For example, for FP4 when n; = 1, no = 1, ki = m1 + Qu, and ke = 
N2 + G2, the matrix A is defined by the equations 


[Ai(1 — Siz, X jn) + Ao(1 — Six.) + w min(Z, m1) + pw min(J, n2)] pi 

= (1 — Xi-1-n,Xn,-1-s)[A1(1 — 80) pi-1,; + WL — 4j2,)min(7 + 1, n2) pij+1] 
+ (1 — Sj0)[A18in,Xng-7 + A2(1 — Xi-1-2,Xm-s)] Dij—-1 
+ (1 — di%,)[(1 — Xi-nyXn.-1-s) Min(t + 1, m1) + N2Xi-n, Fjn.] Di+1,, 


where 


Se lr=m and _jl,r2o 
™ "10, rxm Xr 0, r<0 


and0<i=hk, and0</jS kz. 

We will describe two iterative techniques for solving (3). The first 
method, based on inverse iteration, requires a few expensive iterations. 
The second, based on splitting A into the sum of two matrices, requires 
many inexpensive iterations. 

Inverse iteration may be written as follows: 


668 — THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1981 


Pick p®, an approximate null vector of A. 
Iterate until convergence: 
For k=1,2,---, 


solve Av") = p®-) 


for v™, (4) 
set p® = v/v II. 


Usually only two iterations are required even if the initial guess is a 
random vector. Because of the near singularity of the matrix A repre- 
sented in the computer, the vectors v™ will be very large; they will 
also become richer in the direction of the null space of the A matrix." 

Using a linear equation solver designed for band matrices, eq. (4) 
may be solved in approximately 3k:k2max(ki, ke) locations. When 
kik, > 500, using a sparse matrix solver will require fewer operations 
and less space. However for larger problems (i.e., k: 2 > 1000) it pays 
to use some of the algebraic structure of the problem. For example, for 
FP4 the matrix A and the vectors p and v can be permuted and 
partitioned so that (4) becomes 


BOX\ [v®? pi 
OcY | { v? |} =| pf” 
2QE)] \v)}  \p >] 


where pa corresponds to p;; where 0 <= i < m and 0 </ S ko, pp 
corresponds to p;; where ni <i = ki and n2</j S ke, and pc corresponds 
to pij where 1 = 11,0 S7 S ke andj = nz, m1 <i S ki. The matrices B 
and C have the form 


So Fo Gnj+1 J 
H Si P, K Gn,+2 J 
B= ° . . ’ C= : . ° 2 
AT Sn,-2 Fn,-2 K Gua J 
H Sp,-1 : K Gi, 
where 
p= w(t + )Ta + ixry +, J = —pnily,xq,, 
H = —\iTn41xky+15 K = —hilaxays 
Si = S + yidzsixkot1, and Gi = yilq.xq, + G, 
where 


_JAtpmin(i,m), if i<hi, 
ve | pm, if i=h, 


and S and G have the form 
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Xo €0 Xngtl e 
—)he2 XxX) ei —)2 Xnnt2 e 
S= . . . ; G= ; 
—Az Xky-1 Chy-1 —Az Xky-1 € 
—)2 Xk, —Az Xk, 


where 


xa fhe temin(y, m2), if j< ke, 
i=) uns, if j= he, 


ej = —n min(7 + 1, nz), 
e = —pNo2. 


The matrices X, Y, Z, Q, and E are very sparse (there are only kz + 1 
nonzero elements in X) but they do not possess any relevant algebraic 
structure as B and C do. The eigendecomposition of B can be expressed 
in terms of the eigendecomposition of an m, X m tridiagonal matrix 
and the eigendecomposition of the (k2 + 1) X (Re + 1) tridiagonal 
matrix S. Similarly the eigendecomposition of C can be expressed in 
terms of the eigendecomposition of a gi X q: tridiagonal matrix and 
the g2 X qe tridiagonal matrix G. Using block Gaussian elimination 
and the eigendecompositions of B and C, solving (4) entails formulating 
and solving a dense system of kz + 1 + q; equations. 

The second method, line successive overrelaxation (SOR), is much 
easier to implement and can be easily generalized to more queues. it 
uses the fact that A and p can be permuted and partitioned so that (3) 
can be written as 


To Eo Po 
D, T, Fy Pi 
, = 0, 


Dz-1 Tr,-1 Ex,-1 
Di, = Tr, Pr, 
where p; corresponds to p;;, 0 <j S ke. The T matrices are tridiagonal 
and the E and D matrices are diagonal. The exact elements of the 
matrices depend on the feature package modeled. 


In line sor, initial vectors p,” and a parameter w are chosen and the 
following iteration rule is executed: 


For k = 1, 2, --- until convergence, 
for i = 0, 1,2, +--+, ki, 
solve Tw = Dp}? + Exp, for v, 


set pi"? = pi” + w(-v pM). 
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After convergence, 
set p — p/|p|1. 


The off-diagonal elements are never stored, but generated each time. 
Each iteration requires about 5m multiplications, 7m additions, and 
2m divisions, where m is the size of the state space. The number of 
iterations depends on the ratios of the )’s to p’s, the size of the problem, 
the choice of w, the feature package, and the required accuracy. Most 
of the examples given in Section IV required between 50 and 150 
iterations for problems of between 500 and 2500 unknowns. 

When w = 1, line sor is called block Gauss-Seidel and is known to 
converge (see Kaufman’). For our problems a 10 percent change in w 
sometimes doubles the number of iterations, as Table I illustrates. 
Usually the optimal w is about 1.7 and choosing w above 1.9 causes 
divergence. A heuristic algorithm was developed for adjusting w as the 
iterative scheme progresses. The algorithm uses the theory developed 
in Ref. 12 to obtain an overestimate of w and takes into consideration 
the fact that as w is increased, the relative error tends to initially 
increase before decreasing if the choice of w yields a convergent 
scheme. 

Various other splitting schemes have been tried. Chebychev accel- 
eration with Gauss-Seidel preconditioning requires about the same 
amount of computational effort as line sor with the optimal w.” 
Chebychev acceleration also requires the estimation of certain un- 
known parameters but the course of the algorithm seems to be much 
less sensitive to their values. Chebychev acceleration with an incom- 
plete LU preconditioning has been successful for problems, such as 
the ones we have, in which the graph underlying the matrix is two- 
cyclic. A block sor method in which each diagonal block is a separable 
matrix has also been tried. For the problems we have considered, the 
increase in the amount of work per iteration is not compensated by 
the decrease in the number of iterations. 


Table I—w vs number of 

iterations for FPO, p = 1, 
Ay = 40, Ae = 36, n= 
40, ne = 40, q; = 10, 

G2 = 10, stopping criteria 


~ 107° 
No. of 
ra) Iterations 
1.60 184 
1.70 119 
1.75 94 
1.80 114 
1.90 217 
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IV. NUMERICAL RESULTS 


The criteria we use for engineering the system are that the blocking 
probabilities, L; and Le, and the average delays, W; and Ws, of calls 
which enter the queues, do not exceed specified values. For prescribed 
loads a; = Ai/p and a2 = Ao/p, the aim is to choose the number of 
primary and secondary trunks n; and ne, and the number of waiting 
spaces qi and qo, so that the criteria are met. The procedure used was 
to select values of 7; and no, and then determine the smallest values of 
qi and qe for which the criteria on the blocking probabilities are met. 
The process was repeated for different sets of values of n; and ne, to 
find sets for which the criteria on the average delays are also satisfied. 
A search was made to determine the smallest total number of trunks 
n, + ne for which the criteria on both the blocking probabilities and 
the average delays are met, with appropriate choices of qg: and qo. 

In the numerical results given below, we take » = 1, so that the 
average delays are given in units of mean holding time. Two sets of 
results were obtained. The first set corresponds to primary and 
secondary loads a, = 3 and a2 = 8, and the desired criteria were 
max(Z;, Le) = 0.01 and max(W,, W2) <= 1. The results in Table 
II, which correspond to four primary and nine secondary trunks, 
indicate how changes in the number of waiting spaces affect the 
blocking probabilities and the average delays. The results of the search 
to minimize the total number of trunks required are depicted in Table 
III. In none of the cases were we able to satisfy the criteria with a total 
of only 12 trunks. Moreover, for both Fp8 and FPO, for a total of 13 
trunks we required 4 in the primary group and 9 in the secondary 
group. This was not the case for FP4. However, if we take into account 
the fact that primary trunks are less expensive than secondary ones, 
then for rp4 we would choose the set with four in the primary group 


Table II—Blocking probabilities and average delays for offered loads 
a, = 3 and a = 8, and unit mean holding time, n; = 4 and no = 9 
trunks, and different values of q, and qz2 


FP 71 q2 10°L 1 10°L> WwW, We 
8 6 18 1.063 1.027 0.612 0.727 
8 7 18 0.725 1.027 0.650 0.727 
8 6 19 1.076 0.904 0.613 0.748 
8 7 19 0.735 0.905 0.651 0.748 
0 7 17 1.226 1.070 0.730 0.705 
0 8 17 0.907 1.069 0.777 0.705 — 
0 7 18 1.232 0.943 0.730 0.727 
0 8 18 0.911 0.941 0.777 0.727 
4 6 18 0.945 1.059 0.519 0.727 
4 7 18 0.637 1.060 0.549 0.727 
4 6 19 0.959 0.932 0.521 0.748 
4 7 19 0.648 0.934 0.552 0.748 
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and nine in the secondary group. Also listed in Table III are other 
steady-state quantities of interest, as defined in Section II. In addition, 
Ow = Aili2 + Riz is the total mean overflow rate from the primary 
group to the secondary group. 

It is of interest to compare the results for the three cases in Table 
III corresponding to n; = 4 and nz = 9, bearing in mind the differences 
in the number of waiting spaces and the results in Table I. Calls 
arriving at the primary group are most likely to be queued, and (by 
far) least likely to overflow (immediately), for FPO, as is to be expected, 
since overflow is permitted only when the queue is full. Moreover, the 
average delay in the primary queue is largest for FPO, and larger for 
FP8 than for FP4, since for FP8 no overflow is permitted from the 
primary queue. Also, the total mean overflow rate for rp4 exceeds that 
for FP8, although the immediate overflow probability is smaller for 
FP4. It is apparently the capacity for overflow from the primary queue 
which accounts for the other solutions for FP4 with a total of 13 trunks. 

The second set of results corresponds to primary and secondary 
loads a; = 10 and a2 = 5, and the desired criteria were max(Ii, Le) = 
0.005 and, at first, max(W;i, W2) < 0.8. The results of the search with 
a total of 17 trunks are depicted in Table IV. There are no solutions 
for FPO which satisfy the criteria, and just two for FpP8. Three solutions 
are given for rp4. Although it is not possible to satisfy the criteria with 
a total of 17 trunks and more than 10 in the primary group, it is 
expected that there are solutions with less than 8 in the primary group, 
but we have not checked this. There are no solutions for FP4 or FP8 
with a total of 16 trunks which satisfy the criteria. If we relax the 


Table Ill—Steady-state quantities for offered loads a, = 3 and az = 
8, and unit mean holding time, with max(L,, L2) = 0.01 and 
max(W,, W2) = 1 


FP 8 0 4 4 4 4 

ny 4 4 4 3 2 1 

Ne 9 9 9 10 11 12 

qi 7 8 6 8 10 11 

q2 19 18 19 12 9 7 

10°, 0.735 0.911 0.959 0.955 0.890 0.890 
10°L2 0.905 0.941 0.932 0.831 0.810 0.934 
W, 0.651 0.777 0.521 0.679 0.805 0.861 
W, 0.748 0.727 0.748 0.411 0.284 0.214 
Q1 0.293 0.477 0.288 0.351 0.391 0.413 
Q2 0.682 0.621 0.703 0.563 0.492 0.451 
he 0.080 0.004 0.056 0.142 0.258 0.400 
Ru 0.878 1.430 0.743 0.716 0.545 0.291 
Rie 0 0 0.123 0.336 0.628 0.950 
Ree 5.455 4.967 5.621 4.506 3.938 3.607 
Ow 0.241 0.012 0.289 0.762 1.401 2.151 
Xi 2.737 2.960 2.682 2.209 1.572 0.823 
Xe 8.169 7.937 8.215 8.695 9.336 10.076 
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Table IV—Steady-state quantities for offered loads a, = 10 and 
a> = 5, and unit mean holding time, with max(L,, L2) = 0.005 
and max(W,, W2) = 0.8 


FP 8 8 4 4 4 

nm 10 9 10 9 8 

ne 7 8 7 8 9 

n 20 22 19 20 21 

qe il 8 11 9 7 

10°; 0.481 0.473 0.473 0.480 0.445 
10°L2 0.453 0.481 0.494 0.345 0.430 
WwW, 0.637 0.773 0.508 0.532 0.547 
We 0.460 0.309 0.460 0.317 0.237 
ray 0.433 0.421 0.464 0.474 0.481 
Q: 0.626 0.538 0.683 0.623 0.583 
Ty 0.117 0.187 0.054 0.087 0.123 
Ru 4.327 4.208 3.921 3.622 3.268 
Ri 0 0 0.716 1.120 1.540 
Reo 3.131 2.692 3.417 3.114 2.915 
Ow 1.169 1.875 1.257 1.987 2.774 
x 8.783 8.078 8.696 7.965 7.182 
xX, 6.147 6.851 6.232 6.970 7.752 


criteria on the average delays to max(W;, W2) = 1, then additional 
solutions are obtained, as depicted in Table V. There is now one 
solution for FPO with a total of 17 trunks, and two more solutions for 
FP8. Two additional solutions are given for FP4, one with a total of 
only 16 trunks (and 34 waiting spaces for the primary group). There 
are no solutions for FP8 or FPO with a total of 16 trunks which satisfy 
the relaxed criteria. 

Finally, setting aside the problem of engineering the system, we 
examined the effect on the blocking probabilities and the average 


Table V—Steady-state quantities for offered loads a, = 10 and 
a2 = 5, and unit mean holding time, with max(L,, L2) = 0.005 and 
0.8 < max(W;, W2) = 1 


FP 8 8 0 4 4 

ny 11 8 11 11 10 

ng 6 9 6 6 6 

1 18 24 24 18 34 

Qe 18 7 17 19 19 

10°Z, 0.461 0.482 0.499 0.437 0.466 
10°L2 0.495 0.354 0.456 0.428 0.471 
W, 0.523 0.954 0.753 0.478 0.983 
W, 0.883 0.237 0.866 0.898 0.898 
Q, 0.437 0.406 0.654 0.452 0.676 
Q2 0.761 0.479 0.579 0.795 0.874 
Te 0.054 0.263 0.002 0.025 0.020 
Ru 4.373 4.060 6.540 4,187 6.191 
Riz 0 0 0 0.335 0.568 
Rx 3.807 2.397 2.895 3.973 4.371 
O12 0.540 2.633 0.017 0.587 0.771 
Xi 9.414 7.319 9.933 9.370 9.183 
X2 5.515 7.615 4.994 5.565 5.747 
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Table Vi—Blocking probabilities and average delays for 

offered loads a; and a2 = 0.9a,, unit mean holding time, 

n, = 40 = nz trunks, q; = 10 = qo2 waiting spaces, and 
different values of a, 


FP Qa 10°L, 10°Le2 W, W, 
0 42 5.355 3.654 0.1475 0.1259 
0 40 2.906 2.022 0.1375 0.1162 
0 38 1.291 0.954 0.1270 0.1065 
0 36 0.449 0.378 0.1162 0.0970 
8 42 3.796 4.774 0.1313 0.1259 
8 40 1.758 2.845 0.1184 0.1162 
8 38 0.639 1.436 0.1055 0.1065 
8 36 0.177 0.596 0.0930 0.0970 
4 42 3.588 5.020 0.1179 0.1259 
4 40 1.642 2.995 0.1043 0.1162 
4 38 0.590 1.510 0.0909 0.1065 
4 36 0.161 0.623 0.0782 0.0970 


delays of varying the loads offered to a prescribed system. We did this 
since in a given situation the configuration of the system is fixed, but 
one would expect that the loads would vary in the course of the day. 
We chose a configuration with 40 trunks and 10 waiting spaces in both 
the primary and the secondary groups, to show that our programs can 
handle larger problems than those listed so far. The primary and 
secondary loads were varied simultaneously with a2 = 0.9a1, and the 
results are given in Table VI. As expected, the blocking probabilities 
and the average delays increase with increasing load, the more signifi- 
cant effect being on the blocking probabilities. 

For a given load, the blocking probability for calls from stream S; is 
larger for FPO than for FP8, and slightly larger for Fp8 than for FP4. On 
the contrary, the blocking probability for calls from stream S) is 
smaller for FPO than for FP8, and slightly smaller for FP8 than for FP4. 
These orderings are not surprising in view of the fact that no overflow 
is permitted from the primary queue for either FPO or FP8, and 
moreover overflow is permitted for FPO only when the primary queue 
is full. As before, the average delay in the primary queue is largest for 
FPO, and larger for rP8 than for rp4. The average delay in the secondary 
queue is, of course, the same for FPO, FP4, and FP8. 
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This paper estimates the forecast error distribution for outside 
plant using data from the central office forecast measurement plan. 
We then determine the impact of the forecast errors on feeder-cable 
sizing by using this distribution to estimate the conditional distri- 
bution of engineered cable size with respect to optimum cable size. 
The marginal distribution of optimum cable size is estimated from a 
growth rate distribution which in turn is estimated from cable ship- 
ment data. We then use the resulting joint distribution to weight the 
percentage cost penalty of each possible combination of optimum and 
forecast size. The impact analysis is done separately for each gauge. 
By weighting by the million conductor feet of each gauge shipped, we 
then obtain an estimate of overall sizing-error cost penalty. The 
resulting penalty estimate is about 0.5 percent of the annual feeder- 
cable-construction program. 


l. INTRODUCTION 


A feeder route is a major network of cables extending from the 
central office to within 4 mile or so of customers.’ When a feeder route 
needs relief, a cable size is selected with the goal of minimizing the 
discounted sum of costs over time. Because of forecast errors, however, 
sometimes a cable that is larger or smaller than the optimum is placed. 
One often hears that the feeder-cable sizing curves are so flat that 
sizing decisions are relatively insensitive to forecast errors. On the 
other hand, a small percentage of a large construction program still 
represents a substantial amount of money. In this paper we attempt to 
quantify the impact of forecast deviations on the feeder network by 
first estimating the error distribution and then using it to examine the 
effect of forecast errors on feeder-cable sizing. While it is not presented 
here, a preliminary study indicates that the impact of forecast devia- 
tions on feeder-cable relief timing is at least as great as that on sizing. 
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ll. FORECAST ERROR DISTRIBUTION 


In this section, we derive an estimate for the distribution of forecast 
errors using data from the central office forecast measurement plan. It 
should be noted that all deviations between forecast and actual are 
included here under the forecast error category. Thus the forecast 
errors include some deviations caused by count errors and others 
caused by boundary changes that have not been reflected accurately 
in the records. 


2.1 Nomenclature 

The units of primary concern are available pairs, which include both 
working and idle pairs. The data available for estimating the forecast 
error distribution, however, are in terms of main stations (plus equiv- 
alent main stations). The distribution will therefore be derived first in 
terms of main stations plus equivalent main stations and then con- 
verted into available pairs. 

The basic items of interest are defined here: 


b = base in-service or total value (the actual on which the forecast 
was based), 


t = forecast interval, in years, 
f = forecast in-service value, and 


a = actual in-service value for the date for which the forecast was 
made. 


Several important variables are derived from the above basic ones: 


e = f — a = forecast deviation, 


gr= a ee forecast average of annual growth rate, and 





a-b 
g= = actual average of annual growth rate. 


These items are shown in Fig. 1. 


2.2 Data description 


Several years ago, the central office forecast measurement plan 
(COFMP) was established to collect forecast data from the Bell System 
operating companies. These data were collected for short-term wire 
center forecasts. The data used in this study were collected in the 
fourth quarter of 1978. For each of 1266 wire centers, we had the 
following: 

(t) identification (company, area, and wire center), 
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IN—SERVICE VALUE 





TIME 


€=f-a = FORECAST DEVIATION 
gf =(f-b]/t_ = FORECAST GROWTH RATE 
g =(a—b)/t = ACTUAL GROWTH RATE 


Fig. 1—Definitions of forecast variables. 


(it) b, f, and a, for t = 1, and 

(tit) the number of main stations transferred from a wire center to 
another one during the forecast interval. 

The values b, f, and a are in terms of main stations (plus equivalent 
main stations). 

Table I gives examples of the above data items. 

In addition to the above items, we had two items for each wire 
center that were not used in the study. For each wire center, the 
month of the end of the forecast period was available. Since in all 
cases, the month fell within a three-month period, we did not feel that 
the differences would be significant. A seasonal indicator was also 
available for each wire center. Several wire centers were flagged as 


Table |I—Examples of COFMP data 


Main stations (plus equivalents) 


Com- Wire Se ee eee 
pany Area Center Base Forecast Actual Transfer 
4 17 109 17785 18260 18400 +23 
7 34 344 871 910 891 0 
7 35 368 10277 10964 11021 0 
10 40 430 3171 3240 3296 0 
14 54 584 5715 5925 6013 0 
14 56 596 2187 2365 2286 0 
16 67 1141 40317 42205 42508 +374 
17 68 1158 21618 24626 24378 0 
18 73 1254 2054 2110 2100 0 
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having the annual maximum occur at some time other than the end of 
the year. This information was not needed for the study, since all 
forecast intervals were exactly one year. 

Since it was desirable to model the forecast results in terms of 
growth rates, the forecast and actual values (f and a) were adjusted 
for any wire center that had a transfer by subtracting the signed value 
of the transfer from them. This allowed us to compare the forecast 
and actual growth rates for the original serving areas. 

An initial examination of the data showed that one company had a 
much larger percentage of its wire centers represented than did any 
other company. Of the 1266 wire centers, 507 were from that company. 
The coFMP was intended to collect data only for wire centers that have 
at least 500 main and equivalent-main telephones and that have a 
traffic order prepared while the forecast is in effect. Eliminating wire 
centers that are less than 500 in size, and arbitrarily retaining every 
third wire center of size 500 or greater reduced the representation from 
that company to the point where it was similar to that of the other 
companies. At this point, we retained 863 of the 1266 wire centers, and 
felt that they provided a representative cross section of the Bell 
System. In Fig. 2, the forecast growth rate is plotted versus the actual 
growth rate for these 863 wire centers. 


2.3 Model development 


From Fig. 2, it is obvious that the variance is increasing with g. We 
tried both square root and log transformations and found that the log 
transformation, with a shift of 50 for both g; and g, did an excellent job 
of stabilizing the variance. 

To be able to use a shift of 50, we had to drop 12 data points with g; 
and/or g less than or equal to —50. Six of these points were of relatively 
little interest, for both g; and g were negative, and we expect the model 
to be used for cases where cable is placed to accommodate growth 
(gr > 0) or where it should be placed (g > 0). The other six points had 
either g; or g less than —60, with the other value positive (all except 
one were greater than +120). These six points would have been outliers 
and should have received small weights for any reasonable model. 
Therefore, their loss is not serious. 

The remaining 851 data points are shown: in Fig. 3, where 
In(g; + 50) is plotted against In(g + 50). At a glance, it appears that 
the variance is now decreasing with g. That this is not the case will be 
shown later. Basically, the illusion is due to more points existing at the 
lower g values. 

We used robust regression to estimate the relationship between g; 
and g. After an initial ordinary least-squares step, we used a Huber 
iteration to downweight points with residuals more than 1.5 standard 
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FORECAST GROWTH RATE IN MAIN STATIONS PER YEAR 


0 1000 2000 3000 4000 5000 6000 7000 8000 
ACTUAL GROWTH RATE IN MAIN STATIONS PER YEAR 


Fig. 2—Forecast vs actual growth rates. 


deviations away,’ using the estimate of the standard deviation obtained 
from the first step. We followed the Huber step by a biweight iteration,° 
using a dispersion value six times the Huber estimate of the standard 
deviation. 
The growth rate model is 
In( gj + 50) = a + B In(g + 50) + », (1) 


where a and £ are parameters to be estimated and » is a residual noise 
term with mean 0 and a variance o” to be estimated. The estimates 
resulting from each step of the regression are given in Table II. 
Figures 4 through 7 show residual plots for the residuals from the 
final step. The plots against the dependent and independent variables 
shown in Figs. 4 and 5 indicate reasonably well-behaved residuals. 
Although it would have been difficult to make use of any relationship 
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2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5 9 
In(g + 50) 


Fig. 3—In(g; + 50) vs In(g + 50). 


involving the size of an area for which a forecast is produced, the 
residuals were plotted versus the base size in Fig. 6. Figure 6 indicates 
that there is no structure involving the base size that needs to be 
included in the model. Finally, the residuals are shown for each of the 
19 Bell System operating companies in Fig. 7. Here, too, there is no 
obvious need to include a company effect in the model. The larger 
extreme residuals generally occur for those companies with a larger 
number of data points. 

It should be pointed out that Figs. 4 through 7 show only 847 of the 
851 points. The other four points are shown in Table III. Three of 
these are outliers that received zero weight in the biweight iteration 
and one lies just beyond the range that was plotted. 


Table il—Results of each step of the regression 


Step ri B é R? 
Ordinary Least Squares 0.363 0.929 0.322 0.920 
Huber 0.277 0.943 0.285 0.937 
Biweight 0.250 0.948 0.274 0.942 
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RESIDUAL 


RESIDUAL 





In(g¢ + 50) 


Fig. 4—Residual vs In( g; + 50). 





In(g + 50) 


Fig. 5—Residual vs In(g + 50). 
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Fig. 6—Residual vs In(d). 
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Fig. 7—Residual vs company. 
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MEDIAN GROWTH RATE: In(g + 50) 


Fig. 8—Standard deviation vs In(g + 50). 


To determine if the variance of the residuals is sufficiently constant 
with respect to g, the 851 data points were ranked by g and divided 
into 23 groups of 37 points each. For each group, an unbiased estimate 
of the standard deviation was calculated, using the weights resulting 
from the biweight step in the regression. The resulting values are 
shown in Fig. 8, plotted versus the median values of In(g + 50). No 
overall trend is obvious in Fig. 8 and regression confirms that it is 
reasonable to assume a constant variance. 

A standardized deviate was found for each data point by subtracting 
the value predicted by the regression model (1) and dividing by the 
regression standard deviation, 


_ In(gy + 50) — [0.250 + 0.948 In(g + 50)] 


d 
0.274 


(2) 


A Q-@ plot of these deviates against the standard unit normal showed 
an excellent fit for the bulk of the data but with tails larger than given 


Table Ill—Residuals for the four points not shown in Figs. 4 to 7 
Wire Com- Biweight 
Center Residual In(g;+ 50) In(g + 50) b pany Weight 

128 2.14 4.96 2.71 46,328 4 0 

138 1.86 5.56 3.64 15,232 4 0 

153 —2.17 6.10 8.46 15,587 4 0 

467 1.03 6.52 5.53 17,255 13 0.410 
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(In (g¢+ 50) — [0.250 + 0.948 In(g+ 50}])} / 0.274 





—-4 -3 —2 -1 0 1 2 3 4 
LOGISTIC (alpha = 0, beta = 1.6) 


Fig. 9—Q-Q plot of empirical distribution vs logistic. 


by the normal. Substituting for the normal, a logistic distribution with 
parameters a = 0 and £ = 1.6 (mean 0 and variance 1.29) resulted in 
the satisfactory Q-Q plot shown in Fig. 9. All but three of the 851 
points are shown in Fig. 9. The other three are given in Table IV. 

We conclude that it is reasonable to assume that the density and 
distribution functions of d are given by 


1.6 e164 1 


The modeling was done using coFMP data expressed in terms of 
main stations plus equivalent main stations, so all growth rates used 
so far have been in terms of main stations plus equivalent main stations 
per year. To study the impact of forecast errors on the feeder cable 
network, we need growth rates in terms of available pairs per year. If 
we now define g” to be a growth rate in terms of available pairs per 
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year, and similarly define g) to be the corresponding growth rate in 
terms of main stations plus equivalent main stations per year, the 
relationship can be estimated as follows: 


g” = 0.65 gi? (4) 


The 0.65 ratio in the above expression is the 1977 Bell System average 
main frame fill. Since the actual fill is generally somewhat lower at 
cross sections out on the route than it is at the main frame, this ratio 
tends to slightly understate the effect of forecast errors on available 
pairs. Substituting (4) in (2) for both g; and g and dropping the 
superscripts, gives 


_ In(0.65g; + 50) — [0.250 + 0.948 In(0.65g + 50)] 


d 
0.274 ‘ 


(5) 


where the growth rates are now in terms of available pairs per year. 
From (3) and (5), we find that the conditional density and distribu- 
tion functions of g; with respect to g are 


3.8 eoe4y 
MEN) = 065g; + 50) 4 ey 
and 
1 
F(gl8) =;{ Gr (6) 
where 


y = In(0.65g; + 50) — [0.250 + 0.948 In(0.65g + 50)], 


and where g; and g are given in terms of available pairs per year. 
Expression (6) is used in Section III to estimate the impact of forecast 
errors on feeder-cable sizing. 

First, however, three additional aspects of the model derivation need 
to be discussed. The cormp data are all for a forecast interval of one 
year. How well does (6) work for other values of ¢? One might 
intuitively expect that forecast errors, even when normalized by divid- 


Table IV—Three points omitted 


in Fig. 9 
Wire Empirical Logistic 
Center Value Value 
153 —4,65 —7.92 
138 3.96 6.79 
128 4.65 7.81 
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ing by the forecast interval as is done when dealing with growth rates, 
would be larger for longer forecast intervals. Indeed this tended to be 
the case for 22 wire centers from one company that had data on wire 
centers forecasts for 1-, 2-, 3-, and 4-year intervals. On the other hand, 
similar data from another company shows that the relative error tends 
to decrease for the longer intervals. An experienced forecaster was not 
surprised at this decrease and explained it as follows. It is often easier 
to determine the potential development for an area than to determine 
when an area will achieve that development. In view of the scant and 
conflicting evidence, we decided that (6) should be used for the forecast 
intervals encountered in feeder-cable sizing. 

The forecast error distribution was derived from wire center data 
and is to be used in the feeder-cable network. Wire center forecasts 
are based in part on time series data that often do not exist for portions 
of a feeder route. Thus one would expect that the forecast deviations 
for outside plant forecasts may be somewhat larger than indicated by 
expression (6). In the absence of specific data, however, (6) is used as 
the estimate for outside plant forecast errors. 

Also, if the main frame fill becomes lower than 0.65, the estimated 
forecast deviations in terms of available pairs will be somewhat greater 
than given by (6). 


lll. FORECAST ERROR IMPACT ON FEEDER-CABLE SIZING 


We use the forecast error distribution derived in Section II to 
estimate for each gauge the conditional distribution of discrete-engi- 
neered cable sizes (with size based on forecast) with respect to each 
possible discrete optimum cable size based on actual growth. Important 
assumptions are that growth is linear and that there are no structure 
congestion problems or opportunities to use pair gain systems. If 
inflation is not considered, it would be appropriate to use a 12 percent 
discounting rate. A 6 percent inflation rate for underground cable leads 
to a 6 percent discounting rate when inflation is considered. 

The marginal distribution of optimum cable sizes for each gauge is 
determined from the distribution of growth rate in that gauge, and 
that is, in turn, estimated from cable-shipment data. The conditional 
and marginal cable-size distributions give the joint distribution for 
each gauge of optimum and forecast (i.e., engineered) cable sizes. 

The percentage cost penalty of each possible combination of opti- 
mum and forecast size is determined and multiplied by the probability 
of that combination occurring. These values are summed to give an 
overall sizing-error cost penalty for each gauge. When weighted by 
MCF (million conductor feet) of each gauge shipped, they provide an 
estimate of the cost impact of forecast errors on feeder-cable sizing. 

The following parts of this section describe the steps in detail. 
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3.1 Growth rate distribution 


An estimate of the marginal probability distribution of the optimum 
cable size is needed for each gauge. Were it not for the discounting 
rates used in the past by some companies to size cables, it would be 
possible to estimate these distributions directly from cable shipment 
data. This section estimates the growth rate distribution for each 
gauge, using cable shipment data and an economic sizing relationship 
that relates growth rate to cable size and gauge under the discounting 
rates used. In Section 3.2.1, we use the growth rate distribution derived 
here to estimate, for each gauge, the marginal probability distribution 
of the optimum cable size under the discounting rate that considers 
inflation. . 

1977 shipments of pulp-insulated exchange cable provide the base 
for estimating the growth-rate distributions of feeder cable. The cable 
shipment data give the mcr of each size and gauge shipped. Let F'(x;) 
be the probability that an mcr of pulp cable of the gauge being 
considered is of size less than or equal to x;. 

The economic sizing relationship is used to relate points on these 
cable-size cumulative distribution functions to points on the growth- 
rate cumulative distribution functions. Assuming linear growth, the 
present worth cost of using a cable size, x, to meet a growth rate, g, is* 


b 
PW, g) = SAO (7) 


where 
r = the discounting rate, 
a = the cable intercept cost ($/year/sheath foot), and 
b = the cable incremental cost ($/year/pair foot). 


The values of a and 6 used are the annual charge values for 
underground cables, since most feeder cables are placed in ducts: 


gauge a b 
26 0.38 0.0011 
24 0.36 0.0014 (8) 
22 0.33 0.0020 
19 0.30 0.0036 


It has been estimated that about one third of the companies consid- 
ered inflation for sizing cables that were shipped in 1977. Thus we used 
a weighted average of the two discounting rates to estimate the growth 
rate distribution, 


r= 0.10. (9) 
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Let g(x;) be the growth rate such that the optimum discrete cable 
size is x; for g just less than g(x;) and x; for g just greater than g(x;). 
It can be found by iteratively solving PW(x;, g) = PW(xi+1, g) for g, 
where x;+1 is the next-larger discrete cable size. Substituting from (7), 
&(x;) is the value that yields the following equality: 


(a + bx;:)/r _ (a + bxi01)/r 


Ae ee, 1 
1 — exp[—rxi/g(xi)]_ 1 — exp[—rxis+i/g (xi)] mY 


Because some cable shipments were affected by structure congestion 
and other factors, the raw growth rate distribution found by substitut- 
ing (8) and (9) in (10) is not as smooth as one would expect the actual 
distribution must be. Therefore, a specific form was assumed for the 
distribution, and parameters were determined by fitting the data. We 
found that it is reasonable to assume that[ g(x)]'” is normally distrib- 
uted. Figure 10 shows the data for all four gauges plotted on normal 
probability paper. Using ordinary least-square fits of [g(x)]’” versus 
unit normal standard deviations corresponding to F(x), gives the 
following growth rates: 


Vg ~ N (ng, 63), (11) 
gauge Bee 
26 27.23 8.38 
24 24.08 6.49 
22 19.42 4.77 
19 12.12 3.18 


The above distributions are shown as the straight lines on Fig. 10; 
Fig. 11 shows the density functions for 26, 24, and 22 gauge. As one 
would expect, the growth rates for the finer gauge demand are generally 
larger than for coarser gauge demand. 

One could argue that the cable shipment data more properly lead to 
an estimate of the forecast growth rate distribution, instead of the 
actual growth rate distribution as assumed here. It is not expected that 
this would result in a significant change in the sizing penalty estimate, 
but it could be studied by iteratively assuming distributions for the 
actual growth rate and solving through the conditional cable-size 
distributions until the resulting marginal distribution for the forecast 
cables agreed sufficiently closely with cable-shipment data. 


3.2 Cable-size errors 


For each gauge, the joint distribution of optimum and forecast cable 
sizes is found. Let 


x* = an optimum cable size, 
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Fig. 10—Growth rate distribution. Normal probability plot of the square root of the 
growth rate is based on 1977 Western Electric shipments and 10 percent discounting 


rate. 
x;= a forecast cable size (i.e., an engineered size based on a 
forecast), 
p(x;|x*) = conditional probability of x; given x*, 
Dx*(x*) = marginal probability of x*, and 
p(x*, xs) = joint probability of x* and x;. 
3.2.1 Marginal probability distribution of x* 


The distribution of x* is found from the growth rate distribution of 
Section 3.1. The discounting rate that considers inflation, 6 percent, 
was used to determine the growth rate intervals for which each discrete 


cable size is optimum. That is, 
r = 0.06. (12) 
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Fig. 11—Growth rate density function. Distribution of actual growth rates is based on 
1977 Western Electric shipments and 10 percent discounting rate. 


The growth rate intervals were determined by substituting (8) and (12) 
in (10). The probability of the growth rate being in each of the intervals 
is then found from (11). Table V gives 26-gauge values. Actually the 
interval for the smallest cable size should start at a growth rate of 0, 
but since the probability that g < 0 is so small, it was included with 
that for the smallest cable. 


3.2.2 Conditional probability distribution of x; given x * 


The distribution of x; given x* is found for each gauge from the 
forecast error distribution of Section II. Again, using a 6 percent 


Table V—Twenty-six gauge values 


x g P,-(x*) 
300 —o-— 16.65 0.003 
400 16.65 — 29.97 0.002 
600 29.97 — 61.05 0.005 
900 61.05 — 114.16 0.014 
1200 114.16 — 182.96 0.027 
1500 182.96 — 267.47 0.046 
1800 267.47 — 367.67 0.071 
2100 367.67 — 483.58 0.098 
2400 483.58 — 615.19 0.120 
2700 615.19 — 762.51 0.132 
3000 762.51 — 1007.01 0.186 
3600 1007.01-— © 0.296 
1.000 
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discounting rate to consider inflation, the forecast growth rate intervals 
that correspond to each discrete value of x; being called for are those 
found in Section 3.2.1. This assumes that each cable is sized to 
minimize the present worth cost based on the forecast growth rate. 

The probability that x; is selected when x* is the optimum size is 
then found from the forecast error distribution (6). In that expression, 
g is the value for which x* minimizes (7) when x* is considered as a 
continuous variable. Taking the derivative of (7) with respect to x, 
setting it to zero, and replacing x with x* gives 


* 
ortle_ TH" _ ar 


g bg” 


A value for g is found for each gauge and x* such that the above 
equality holds when (8) and (12) are substituted for a, b, and r. For a 
26-gauge example, let x* = 1200. The corresponding growth rate, g, is 
148.6 available pairs per year. Figure 12 shows the distribution of gy, 
given this value of g, and the intervals corresponding to each cable 
size being selected. Figure 13 then shows the distribution of x;, given 
x* = 1200. 

Since one would expect the cost penalty for not placing a cable when 
one should have been placed to be at least as large as the penalty due 
to placing the smallest cable available, the probability that g; is 
negative is added to that for the smallest x;. 


3.2.3 Joint probability distribution of x* and x; 


The joint probability is the product of the marginal probability of 
Section 3.2.1 and the conditional probability of Section 3.2.2. That is, 


P(x*, Xf) = px-(x*) p(xp|x*). (13) 





g¢ (N AVAILABLE PAIRS PER YEAR 


Fig. 12—Conditional density of g; given g. 
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Fig. 13—-Conditional density of x; given x*. 


3.3 Penalties for size errors 
Let C(x*, x-) denote the percentage cost penalty when x; is selected 
and x* is optimum. Using the present worth notation of (7), 
PW(x;, 8) — PW(x"*, g) 
PW(x*, g) 
As in Section 3.2.2, g = g(x*) is the value for which x* minimizes 


PW(x*, g), for each discrete value of x*. Expression (14) is evaluated 
by substituting (7) for PW(x, g), (8), and (12) for a, b, and r. 


C(x*, x) = 100 (14) 


3.4 Overall cost of size errors 


For each gauge, the overall penalty is the sum of the cost penalties 
(14), weighted by the probabilities (13). That is, the expected penalty 
for one gauge is 


YL p(x", x)C(x*, x7). 


x* x 


The overall penalty is the sum of the above penalties weighted by the 
1977 Western Electric shipments of pulp-insulated exchange cable of 
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Table Vi—Penalty for each gauge and overall weighted 


penalty 
Percent Siz- Contribution 
ing Error to Overall 
Cost Pen- Cable Ship- Sizing Error 
Gauge alty ments (%) Cost Penalty 
26 0.534 44.6 0.238 
24 0.462 37.0 0.171 
22 0.502 18.2 0.091 
19 0.542 0.2 0.001 
100.0 0.501% Expected overall 


cost penalty 


each gauge. Table VI gives the penalty for each gauge and the overall 
weighted penalty. 

Figure 14 shows graphically the main factors that contribute to the 
sizing-error cost penalty of 26-gauge feeder cable. The solid curves give 
contours of equal cost penalty, C(x*, x;). The dashed curves give the 
1, 10, 90, and 99 percent points on the cumulative conditional distri- 
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Fig. 14—Factors contributing to sizing-error cost penalty of 26-gauge feeder cable (6- 
percent discounting rate). 
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bution of p (x;|x*). Vertical lines show the 1, 10, and 90 percent points 
on the cumulative mcr distribution of cable shipments of 26-gauge 
pulp cable. The cross hatches emphasize the area of greatest interest, 
where 80 percent of the shipments occur and there is an 80 percent 
chance of finding x;, given x*. The cost penalty is less than 1 percent 
in most of this region. 


IV. SUMMARY 


We have derived an estimate for the outside plant forecast error 
distribution. We give this distribution, expression (6), as the condi- 
tional distribution of the forecast growth rate with respect to the 
actual growth rate. We then used it with cable-shipment data to 
estimate that the feeder-cable sizing penalty due to forecast errors is 
about 0.5 percent of the annual feeder-cable-construction program. 
This represents a substantial amount of money, even though it is a 
small percentage, because of the large construction program. 
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The Effects of Misclassification Error on the 
Estimation of Several Population Proportions 


By J. D. HEALY 
(Manuscript received October 21, 1980) 


Assume that each item from a set of items 1s classified into one of 
several categories. We then use the proportion of items classified into 
a category to estimate the true proportion of items from the category 
in the population. This article models the effect of misclassification 
error on the estimate of the true proportion. We discuss two conditions 
which can be used to determine the adequacy of a classifier. We 
present an optimal classification algorithm which can be used when 
the joint distribution of the variables on which classifications are 
based is known separately for items from each category. 


1. INTRODUCTION 


Let us suppose that we observe a set of items which can be split into 
several distinct categories. Each item is measured and classified by 
some device into one of these various categories. However, the classi- 
fied category for an item and the true category may not be the same, 
ie., the device may make a misclassification error. The observed 
proportion of items in a category is then used to estimate the true 
proportion. 

The preceding scenario often occurs in quality control’ and medical 
research.”* In quality control, individual manufactured items from a 
sample or lot are often classified by a mechanical device as defective 
or not and the proportion of defectives in the sample is then used to 
estimate the proportion of defectives from the entire process. In 
medical research, the items are people and the idea is to estimate the 
proportion of people with various diseases. In a Bell system example, 
the items would be phone calls and the categories would be busies, 
completed calls, reorders, etc. An automated device would attempt to 
determine the true category for each call. The output of the device 
would then be the estimated proportion of calls in each category. This 


697 


last example motivated the analysis contained in this article. Note 
that the problems discussed here are quite different from the tradi- 
tional classification problem,’ where the goal is to maximize the 
probability of correctly classifying each item. In all the above examples, 
errors from misclassification can have a serious effect. 

In this article, we attack two separate problems. In the first problem, 
we assume we have little or no control over the internal design of the 
classifier; all that we have is an estimate of the probability of classifying 
items from each category into each of the other categories. The object 
is to develop a simple way to specify how good the classifier must be. 
Also, we should indicate to the designer the direction in which im- 
provements are necessary. The second problem handles the case when 
we do have control over the design of the classifier. In this case, we 
assume that object is to design the classifier so that the effects of 
misclassification error are minimized. In this article, we are concerned 
only with a classifier’s ability to estimate proportions, i.e., our loss 
function is entirely different than the usual loss function. 

This article is organized in the following way. Section II introduces 
notation and explains the effects of misclassification error on the 
estimated proportion of items in a category. Section III discusses the 
case when we have little or no control over the design of the classifier. 
In Section IV, we discuss the case when we have control over the 
classifier. The resulting minimization problem involves a function that 
is quadratic in the probabilities of classifying items. 


ll. EFFECTS OF MISCLASSIFICATION 


Let m be the number of categories. All vectors in this paper are m- 
dimensional column vectors, and matrices are m by m matrices. Also 
let p, p*, and p be the m-dimensional vectors of the true probabilities 
of items from the various categories, the probabilities of classifying 
items into various categories, and the observed proportion of items 
actually classified into different categories, respectively. Let A = (a,j) 
be the m by m misclassification matrix which contains the conditional 
probabilities of classifying items into different categories. For example, 
ay. would be the probability of classifying an item from category 2 as 
an item from category 1. Each column of A sums to one. The diagonal 
elements of A are the probabilities of correct classification, and the 
off-diagonal elements give the probabilities of misclassification. A 
perfect classifier would have A = J, the identity matrix. By the law of 
total probability, p* = Ap. 

In measuring the effectiveness of p as an estimator of p, we use the 
matrix of mean-squared errors. The matrix of mean-squared errors 
contains the mean-squared error of the individual terms and also the 
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cross product terms which indicate how the different errors are related. 
The m by m matrix of mean-squared errors (MMSE) is 


E[(p — p) (p — p)’|p] 
= cov(p|p) + [E(p|p) — pJ[E (|p) — p]’,  () 


where “cov” means the covariance matrix. The diagonal of the covar- 
iance matrix measures precision of an estimator while the diagonal of 
the second term in (1) is the bias squared. 

We assume that 7 items are to be classified. The expected number 
of items from each category is np. If items are classified independently, 
the distribution of np will be a multinomial distribution with param- 
eters Ap and sample size n. From Ref. 6, we obtain the cov(p|p), 


cov(p|p) = [D* — App’A’//n, (2) 


where D* is a diagonal matrix with diagonal equal to p*. The cov(p|p) 
can be separated into two parts, one part which is the covariance 
matrix if the classifier were perfect, and the second part which is an 
adjustment in the covariance matrix because the classifier is not 
perfect: 


cov(p|p) = [D — pp’]/n + [D* — D— App’A’ + pp’]/n, 


where D is a diagonal matrix with diagonal equal to p. 
Since E'(p|p) = Ap, the bias term in (1) becomes 


[E(p|p) — p][E(b|p) — pl’ = (A — D)pp’(A — 1)’. (3) 


Note that this term is not divided by n. Putting the above statements 
together, (1) becomes 


E[(p — p)(p — p)’|p] = (D — pp’)/n 
+ (D* — D— App’A’ + pp’)/n + [((A —I)pp'(A— 1’). (4) 


This equation includes the sampling-error effect (first term on right) 
and the effect of a misclassification error (other terms on right). 


ll. SPECIFYING ACCURACY WITH LITTLE CONTROL OVER THE 
CLASSIFIER 


Assume that our information on a classifier is confined to the matrix 
A, i.e., someone else is responsible for the design of the classifier. We 
may influence the form of A but we have no control over the actual 
functioning of the classifier. Suppose we have a good idea of how large 
the mean-squared errors for each category estimate [the diagonal 
elements of (4)] can be for the application of the classifier. We need 
several guidelines on the form of A which will insure that the classifier 
and its resulting mean-squared errors are adequate. Clearly, we do not 
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want to put constraints on every element of A. Also, we cannot tell the 
designer of the classifier that the mean-squared errors of his classifier 
are inadequate without providing some guidance on how to improve 
the classifier. 

Intuitively, we would like to pick A so that the bias term in (4) 
disappears. We cannot pick A so that the bias term in (4) disappears 
for every p. For each A, however, there is some value of p for which 
the bias term disappears; in fact, A tends to produce p, which are 
collapsed toward the value of p for which the bias term disappears. A 
reasonable strategy is to pick A so that the bias term disappears for 
our “best guess” for p which we denote by po. This means A should be 
picked so that (A — J)po = 0, i.e., po is approximately an eigenvector 
of A with eigenvalue 1. There are many A which satisfy (A — I)po ~ 
0, but which have large mean-squared errors for values of p near po. 
To ensure that mean-squared errors are small for values of p near po, 
we must additionally require that the a;; (diagonal elements of A) be 
reasonably large; note that if the eigenvector condition is nearly 
satisfied, the requirement on the a,; may, in many cases, be quite loose. 

These two conditions: (1) (A — I)po = 0, and (2) ai large, are 
generally easy to check. Assume we have a set of s items which we 
know contain spo items from the respective categories. These items 
are classified and are the results used to estimate A. The first condition 
states that the number of items classified into a category is roughly 
equal to spo. If this condition is not originally satisfied, the designer 
can usually satisfy it by adjusting several thresholds which determine 
where. the classifier places items. The second condition says that the 
classifier cannot misclassify a high proportion of items from any one 
category. The designer is then told which categories do not satisfy this 
condition. . 

We now show that these two conditions can be justified analytically 
when we assume that p has an underlying Dirichlet distribution. That 
is, we now allow p to vary, for example, with environment. The 
Dirichlet distribution is the natural multivariate generalization of the 
beta distribution and it is the conjugate prior for the multinomial 
distribution. We pick the parameters of the Dirichlet distribution so 
that 


E(p) = po, (5) 
cov(p) = —popo/(v + 1) + Do/(u + 1), (6) 


where Dp is a diagonal matrix whose ith diagonal element is the ith 
element of po, and uv is a parameter that indicates how spread out the 
Dirichlet distribution is. These parameters, po and vu, can be chosen so 
that the resulting Dirichlet distribution models the expected environ- 
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mental variability of p. If we take expected values of the terms in (4) 
using (5) and (6), we obtain 


E(p — p)(p — p)’ = Di/n — ADA’/n(v + 1) 
— ApopoA’v/n(v + 1) + (A — I)popo(A — I)’v/(u + 1) 
+ (A —I) D(A —I)’/(v+ 1), (7) 


where D3 is a diagonal matrix whose diagonal equals Apo. This 
equation incorporates the effects of varying p. If n and vu are at all 
large, the fourth term in (7) will be important. If (A — I)po = 0 holds, 
this term will drop out. The fifth term is minimized if the a;; are large. 
Since the second and third terms are generally unimportant [they are 
divided by n(v + 1) which should be large] and since the first term is 
present even with a perfect classifier, satisfying the two conditions will 
minimize the effects of misclassification. In short, we have presented 
a way to require the A matrix to be “near” the identity matrix without 
putting constraints on each and every element of A. 


IV. DESIGNING A CLASSIFIER THAT MINIMIZES MEAN-SQUARED 
ERROR 

Assume that our job is to develop a classifier that minimizes the 
effects of mean-squared error. More specifically, assume we measure 
a vector of variables (x) on each item. Regions R; are defined such 
that if x € R;, we classify the item into category 1. We want to define 
the R; that minimizes a weighted sum of the mean-squared errors for 
the various categories, i.e., that minimizes 


tr[Q E(p — p) (p — p)’), (8) 
where Q is a known positive-diagonal matrix, tr is the trace operator, 
and E(p — p)(p — p)’ is defined by (7). The matrix Q is just a 
weighting factor that can be used to emphasize the important cate- 
gories. Minimizing (8) is equivalent to minimizing 

tr(QABA’ — AC), (9) 


where 
_ a 1 
‘ n(v + 1) 


2 1 
C= (sy (Do + Popov) — = po(1, 1, «+. je. (11) 


(Do + Popov), | (10) 


We handle a more general case by allowing B to be any known 

symmetric, nonnegative definite matrix, and C any known matrix. 
Equation (9) is interesting because it is quadratic in A. The usual 

Bayes multiple decision rule’ is a special case of (9), since it is the rule 
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that minimizes (9) when B = 0, and C is a diagonal matrix with 
diagonal equal to the vector of prior probabilities. Classical discrimi- 
nant analysis’ is a special case of the Bayes multiple decision rule for 
which the distribution of x when the item comes from category 2 is 
multivariate normal with mean p; and covariance matrix 2. 

Let f;(x) be the density of x if x is measured on an item from 
category i. The optimal classification algorithm is given in the following 
theorem: 

Theorem 1: Equation (9) will be minimized if x is classified into 
the ith category when the ith element of 


(fi(X), f2(X), «++ , fe (%)) (2BAQ’ — C) (12) 


is the smallest. 
(Proof: See appendix.) 

In general, applying Theorem 1 should be quite difficult since A has 
to be solved for in (12). We now discuss the two-category case and 
then specialize to the case when x has a multivariate normal distri- 
bution. We give a simple iterative procedure to calculate the required 
quantities for this last case. 

For the two-category case, Theorem 1 reduces to the following 
corollary. 

Corollary 1: If there are only two categories of measurements, then 
eq. (9) will be minimized if x is placed into category 1 if 


fi(x)/fa(x) > K, (13) 
where 
K = [220 (qui + G22) Ger — 2621 (Qui + G22) Ai2 
+ 2(b21G@11 — b22GQ22) + (C22 — C21) ] (14) 
/[2611 (qu + G22) Qi2 — 2b12(Qui + G22) dar 
+ 2(b12G22 — bugu) + (e1 — ¢r2)], 


and bj, qi, and cj are elements of B, Q, and C, respectively. 

Let us now assume that, if x is an observation from category 
i (t = 1 or 2), it has a multivariate normal distribution with known 
mean vector p; and known covariance matrix 2. Then (13) becomes 


log (fi (X)/f2(X)) = X'Z7* (a — pha) — 1/2 (pr 
+ f2)"Z"" (my — fz) > log K. 


As in classical discriminant analysis,’ the distribution of log (fi (x)/ 
f2(x)) is normal with mean a/2 or — a/2 when x comes from categories 
1 or 2, respectively, where 
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= (jy — pe)/Z7* (wa — pe). (15) 


In either case, the variance of log (fi (x)/fe(x)) is a. Therefore, when x 
comes from category 1, 


[log (fi(x)/fa(x)) — 0/2)/Va 
has a standard normal distribution. Similarly, 
[log (fi(x)/fa(x)) + a/2]/Va 


has a standard normal distribution when x comes from category 2. The 
misclassification probabilities can be defined in terms of a standard 
normal random variable, 


ax, = P (z <te) (16) 
ay2=P (z > we Tel) ; (17) 


where Z has a standard normal distribution. Using (14), (16), and (17) 
the values of @i2, doi, and K may be calculated iteratively: 

(1) Let K = 1 and calculate aj2 and a2; using (16) and (17); 

(2) Obtain a new value of K by substituting aj. and az: into K; and 

(3) Repeat the entire process. 

In summary, assume we are trying to minimize tr[QE(p — p) 
(p — p)’], where @ is a known weighting matrix. Also assume we have 
only two categories and that if x comes from category i, it has a 
multivariate normal distribution with mean vector yp; and covariance 
matrix &. The parameter a may be calculated using (15). The B and C 
matrices should be calculated using (10) and (11). Equation (18) gives 
the decision rule for classification into one of the two categories, where 
the values of K, ai2, and a2; are obtained in an iterative manner using 
(14), (16), and (17). 

To apply any of the preceding theory, some prior knowledge of the 
distribution of p is required to estimate po and v. If the parameters jn, 
fiz, and & are unknown, they may be estimated from a sample of data 
with the usual sample means and pooled covariance matrix. An esti- 
mator of the parameter a could then be calculated using (15) with the 
estimates of 1, fl2, and & substituted into (15). The algorithm discussed 
above could then be used to generate A and K where the estimator of 
a is used in (16) and (17) instead of a. The properties of the procedure 
when estimators of the parameters are used require further study. 


V. SUMMARY 


This article presents a model that incorporates the effects of mis- 
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classification error. Several guidelines are presented which can be used 
to determine if a given classifier is adequate. For the case when the 
classifier is yet to be designed, we have given an optimal classification 
algorithm. We discuss the two-category case in detail. 


APPENDIX 
Proof of Theorem 1 


Let R? be the regions that result if the theorem is applied. Let Ao be 
the resulting misclassification matrix. Let R} and A; be the correspond- 
ing elements for some other decision rule. Now consider (9) evaluated 
at A, minus (9) evaluated at Ao: 


tr (QAi BA‘ — QAoBAo — AiC + AoC) 
= tr(A; = Ao) B(Ai = Ao)’Q + 2tr(Ai = Ao) BAoQ 
— tr(A; — Ao)C. (18) 


The first term on the right of (18) is nonnegative since B and Q are 
nonnegative definite. We still have to show that the rest of (18) is 
nonnegative. Consider 


2tr(Ai — Ao) BAoQ — tr(Ai — Ao)C = trAiE — trAok, (19) 


where F = 2BA0Q — C. Let e;; be the i, jth element of E. Equation 
(19) now becomes 


tr(A,F) — tr(AoF) = » $1 (i| x) f(x) adx-e; 
= 2 | eotht20f (X) AX+ emp 
= [ eoecein 


: b f(x) ej: — Y f(x) ent | dx, (20) 


where 
Sinn | Le xER}, 
oe c xER?, 
+e, Jl xERi, 
se alae 1 xER}, 
Since ¢o (k| x) will be zero whenever [)); fj(xeji — Yim fm€ mr] is negative, 


(20) is always nonnegative. Q.E.D. 
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Adaptive Post-Filtering of ADPCM Speech 


By N. S. JAYANT 


(Manuscript received December 3, 1980) 


We show that the quality of adaptive differential PcmM (ADPCM) 
speech can be significantly improved by passing it through a recon- 
struction low-pass filter that is matched to an appropriately defined 
short-time speech cutoff frequency. Practically, the adaptive proce- 
dure involves switching the decoder output into one of a bank of N 
low-pass filters whose cutoff frequencies span the expected range of 
input speech bandwidth. For the case of equally spaced filter cutoffs, 
and with uniform probability density function models for the quan- 
tization noise spectrum and the cutoff frequency, more than one-half 
of the maximum adaptive filtering gain is realizable by a bank of 
four filters. Computer simulations of 16- and 24-kilobit/s ADPCM 
coders indicate that perceived quality gains are in fact greater than 
what is indicated by an analytically predicted objective gain of 2.6 
aB. 


1. SHORT-TIME CUTOFF FREQUENCY 


ADPCM (adaptive-quantization/differential PcM) speech coding usu- 
ally assumes a time-invariant model for speech bandwidth (such as 
3200 Hz for telephone quality applications), and a corresponding time- 
invariant low-pass filter with cutoff fo = 3200 Hz for the decoded 
speech. However, short-time speech spectra of 3200-Hz-filtered speech 
exhibit cutoff frequencies f,, that are often significantly smaller than 
the long-time nominal cutoff frequency fo.’ Figure 1 sketches a short- 
time spectrum at time ¢, and defines f.(t) = f7(t) as the low-pass cut- 
off frequency that includes all but T percent of short-time spectral 
energy. Figure 2 shows long-time-averaged spectra for four sentence- 
length 3200 Hz band-limited utterances [L, B, C, and D], denoting [“A 
lathe is a big tool” (female utterance), “An icy wind raked the beach” 
(female utterance), “The chairman cast three votes” (female utter- 
ance), and “This is a computer test of a digital speech coder” (male 
utterance) ], respectively. It is clear that all the four spectra are low 


707 
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0 felt)= #7 (t) 


Fig. 1—Definition of short-time cutoff frequency f(t). The shaded area includes T 
percent of short-time speech power. 


pass in a long-time-average sense, a fact that is well exploited in fixed- 
prediction differential coding.” Figure 3 shows corresponding histo- 
grams of short-time cutoff frequency fZ(t) for a threshold T = 1 
percent. It is seen that on a short-time basis, speech segments can be 
either low pass (say, fx(f) < fo/2) or all pass (say, fx(t) > fo/2), 
although both of these segment types come from inputs that are low 
pass from a long-time-averaged energy viewpoint. It is also clear that 
the four histograms of Figure 3 are very different; however, as a single 
descriptor of these histograms, we propose a uniform probability 
density model for fz(t), 


. PUT) =F, O< full) < fo (1) 
The adequacy of the above model is clearly a function of the threshold 
T. Clearly for the extreme cases of JT'= 100 percent and T'= 0 percent, 
the pdf of f2Z(t) would degenerate into delta functions at f = 0 and fo, 
respectively, with corresponding low-pass counts of 100 and 0 percent. 
The uniform density model (1), on the other hand, implies equal, 50 
percent occurrences of low-pass and all-pass segments, and Fig. 4 
shows that this is reasonable as a sentence-ensemble average for 
thresholds T = 1 and 2 percent. A threshold of JT’ = 1 percent has been 
used in all of our ADPCM simulations. 

The value of T = 1 percent produces a filtering distortion of 
10 log(1/100) = —20 dB. More relevantly, it constitutes a “perceptually 
acceptable” low-pass threshold, with a filtering distortion that is ob- 
vious only in critical listening. On the other hand, at T = 2 percent, 
the filtering distortion begins to get obvious, and undesirable even for 
low-quality application such as 16 kilobit/s. The role of threshold T is 
discussed at greater length in a recent work which relates to the 
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refinement of ADM (adaptive delta modulation) with adaptive post- 
filtering." This work also shows that consideration of the low-pass 
threshold f7(t) is more useful, for purposes of reducing coder noise, 
than the consideration of a high-pass threshold f #(¢); in other words, 
providing a band-pass reconstruction filter matched to short-time 
high-pass and low-pass cutoffs [fx(t), fxr(t)] was not significantly 
better than providing a low-pass reconstruction filter matched to the 
range [0, fx(t)]. 
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_ Fig. 2—Long-time averaged spectra for four sentence-length inputs. A 3.2-kHz band 
limitation is obvious in all four examples. 
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_ Fig. 3—Histograms of short-time cut-off frequency Z(t) (T = 1 percent) for the four 
inputs of Fig. 2. The short-time cutoff is in general significantly less than 3.2 kHz, the 
cutoff in the long-time-averaged spectra of Fig. 2. 


ll. ADAPTIVE LOW-PASS FILTERS 


Figure 5a shows the effect of a low-pass reconstruction filter in ideal 
adaptive post-filtering; out-of-band ADPCM quantization noise compo- 
nents in the cross-hatched range [fz(t), fo] in the noise spectrum are 
rejected by an ideal low-pass filter matched to f(t). 

Figures 5b and 5c show suboptimal but practical versions of 5a, 
where the value of f(t) causes switching the decoder output into one 
of a bank of N low-pass filters, leading to a degree of noise-rejection 
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Pr [FZ1n< 1600 Hz] PERCENT 





0.5 1 2 
T IN PERCENT 


Fig. 4—Percentages of segments with f7Z(t) < fo/2 in the four inputs of Fig. 2. For 
T = 1 and 2 percent, the mean percentage of such segments is in the order of 50 percent. 


that is always less than that in Fig. 5a. For example, in the N = 2 
example of Fig. 5b, the filter bank consists of two filters with fixed 
cutoffs f.2(= fo) and f.1; in the upper illustration in Fig. 5b, the value 
of f1(t) is not small enough to switch in the lower filter f.1; conse- 
quently there is no out-of-band noise rejection similar to that in the 
upper example of Fig. 5a. With the uniform pdf model (1), the two- 
filter bank system realizes out-of-band noise rejection only 50 percent 
of the time (when f1(t) < fa = fo/2). The four-filter system of Fig. 5c 
is clearly more effective; it realizes nonzero noise rejections in three 
out of the four cases shown, and indeed for 75 percent of all speech 
segments, if the uniform pdf (1) is valid. 

A block diagram of an N-filter-bank adaptive system appears in Fig. 
6. Note that for simplicity, the cutoff frequencies are equally spaced, 


fen = for, (2) 


and that filter n is switched on when the input frequency cutoff is in 
the appropriate ( fo /N)-wide range, 


n-1e fin (3) 


Switch to filtern if N ra N 








Ill. SEGMENTAL s/n GAINS G(N) 


The spectrum of quantization noise depends on many factors, in- 
cluding the nature of adaptive quantization, the input spectrum, and 
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Fig. 6—Block diagram of adaptive low-pass filtering system with an N-filter bank. 
Allowed cutoff frequencies are equally spaced, being integral multiples of fo/N. The 
entire dashed box constitutes an optional refinement of a conventional Appcm decoder. 


the effect of predictor; but experience indicates that the white-noise 
spectrum of Fig. 5 is indeed the single most reasonable model.’* 
Combining this assumption with the uniform pdf for the cutoff fre- 
quency f7(t), we shall now develop expressions for expected gains in 
segmental” s/n due to adaptive post-filtering. 

When the input cutoff f:(t) is such as to switch filter n (n = 
1, 2, ---, N), the noise rejection factor is N/n (see Fig. 5), with a 
maximum value of N for the extremely low values of fx(t), and a 
minimum value of 1 for extremely high values of f,(t). The expected 
gain (in dB) is therefore given by 


N 
G(N) = 3} (10 log *) . (Pr [n] = x) dB, (4) 
n=] 


where the probability of switching in filter n is Pr [x] = 1/N, a result 
of the uniform pdf (1). The filtering gain G(N) can be simply rewritten 
in the form 


10 ¥ 
G(N) = 10 log N —W ¥' log n. (5) 
1 


Figure 7 plots G(N) as a function of N. 
The asymptotic value 





f(t) 
can be evaluated simply by using the identity 


fo 
G(o) = | (10 tog fo )- pUfetnafut (6) 
0 
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Fig. 7—Segmental s/n gain G(N) versus N. This characteristic assumes uniform pdf 
models for coding noise and f(t). The asymptotic value is G(o) = 4.35 dB. More than 
half this gain is realized with N = 4. 


fiInz=zlnz-—z, 
and this results in 
G(«) = 10/e = 4.35 dB. (7) 


This asymptotic value is indeed close to the ideal adaptive filtering 
gains reported in earlier-cited ADM experiments at several bit rates.’ It 
is also seen from Fig. 6 that the four-filter bank method of Fig. 5c 
theoretically realizes more than one-half of the maximum possible dB 
gain G(o) in segmental s/n: G(4) = 2.6 dB. 

The above analytical formulation can also be used to assess the 
efficiency of the equispaced filter design in (2) and Fig. 6. For illustra- 
tion, G(2) = 1.5 dB with N = 2. It-can be shown analytically that an 
optimal design, one that maximizes adaptive filtering gain for N = 2, 
is one for which f.. = (fo/e), rather than f)/2. Figure 8 plots the 
theoretical expected gain in segmental s/n as a function of fa. The 
maximum gain is 1.6 dB, only 0.1 dB better than that in the simple 
design of (2), which suggests fe = fo/2 for N = 2. 


IV. RESULTS OF COMPUTER SIMULATIONS 


ADPCM coders with two- and four-bank filtering systems were simu- 
lated with the speech inputs mentioned in Section J. The coder bit 
rates were 24 and 16 kilobit/s, corresponding to 3- and 2-bit/sample 
quantizers, and an 8-kHz sampling of the inputs. The quantizers were 
adaptive. Both backward adaptive (AQB) and forward-adaptive (AQF) 
quantizers were simulated;* filtering gains were evident in both cases, 
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Fig. 8—Segmental s/n gain wit ‘N= 2 as a function of fcz/ fo. The maximum value 
of 1.6 dB occurs for fe/fo = e~’. The scheme of Fig. 6 suggests fa/ fo = 0.5 for 
N = 2. This produces an s/n gain Realy 0.1 dB less than the maximum value of 1.6 dB. 


but in terms of absolute quality, the A@F coders were clearly better, 
especially at 16 kilobit/s. AQF coders, however, require the explicit 
transmission of step-size information; using for example, four bits once 
for every segment of 256 samples (16 ms). 

The cutoff frequency f1(t) was computed once for every segment of 
256 samples. These segments were Hamming-windowed, zero-padded 
for better frequency resolution, and analyzed by means of a 512-sample 
FFT. 

The filter banks consisted of 33-point Fir filters whose frequency 
responses are shown in Figure 9. Although gains due to a two-filter 
system were almost always noticeable, they were not always signifi- 
cant; specifically they were not significant for inputs with predominant 
all-pass (> fo/2) segments; see Figs. 3 and 4. A four-filter bank, on the 
other hand, provided significant quality improvement with all the four 
test inputs. The perceived improvement in quality was much greater 
in all cases than what the theoretically predicted objective gain of 
G(4) = 2.6 dB indicates. The measured gains were very input-depend- 
ent with an average value slightly less than the theoretical value of 
G(4). This result reinforces earlier results for time-varying filtering of 
ADM speech,’ where again perceptual gains were in excess of objectively 
measured segment s/n improvements; this is probably because the 
residual in-band noise after adaptive filtering is much less annoying 
because of masking by the input signal. Adaptive bandwidth ADPcM, 
of course, has the additional possibility of variable bit allocation. For 
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Fig. 9—Frequency responses of low-pass filters in an N = 4 filter bank. Each filter is 
a 33-point FiR design, and the cutoffs correspond to those in Fig. 5c and Fig. 6, with 
f o = 3.2 kHz. 


example, significant quality improvements have been noticed in a 
system where the input was subsampled at 5.33 kHz whenever 
f(t) < 2.4 kHz, and more quantization bits were allocated for subsam- 
pled segments. For example, in 16 kilobit/s operation these segments 
were coded with three-bits/sample instead of two. This type of variable 
bit allocation is ruled out by definition in aDM which is a one-bit/ 
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sample system. Notice finally that appcm with out-of-band noise 
rejection and variable bit allocation is very similar to frequency-do- 
main subband coding.” 


V. CONCLUSIONS 


We have indicated a conceptually simple but very effective proce- 
dure for improving the quality of ADPCM speech coding. The extra 
information for effecting this improvement is very little; with a four- 
filter design this extra information would be logs 4 = two bits, corre- 
sponding to four possible ranges of cutoff frequency f(t). The com- 
plexity involved in terms of computing fz(t) is quite significant in 
relation to the simplicity of the conventional, basic appcm coder. The 
increased complexity, however, will be relatively less objectionable in 
voice storage applications than in transmission systems; and in either 
case, an attractive feature is that the post-filtering procedure (the 
boxed portion of Fig. 6) can be used as an entirely optional refinement. 
Also, as noted earlier, the adaptive post-filtering technique discussed 
in this paper can be used with significant gains in the context of coders 
other than ADPcM.' 
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Variable Bandwidth Adaptive Delta Modulation 


By J. O. SMITH* and J. B. ALLEN 


(Manuscript received December 3, 1980) 


The ADM(adaptive delta modulation) speech coder is generally used 
with a time-invariant low-pass filter at the decoder output. The 
purpose of this low-pass filter is to reject coder noise at frequencies 
above the. fixed speech passband. The speech spectrum, however, 
tends to occupy only the lower frequencies within the passband during 
voiced speech, and is somewhat “high pass” during unvoiced speech. 
In this paper, we show how the quality of ADM may be significantly 
improved by adaptively filtering the coder output such as to follow 
the natural bandwidth of the speech. This was found to reduce 
drastically the perceived noise in the output of the ADM coding system 
at low bit rates. The use of an adaptive low-pass filter realizes almost 
all of this quality gain. (An adaptive high-pass filter seems to reject 
less audible noise components and seems more prone to introducing 
objectionable artifacts.) We also discuss a method for reducing the 
bit rate with little or no sacrifice in quality (relative to normal ADM) 
by adapting the sampling rate along with the time-varying low-pass 
filter. . 


1. INTRODUCTION 


In this paper, we explore two methods for better utilizing the time- 
varying bandwidth of speech in ADM (adaptive delta modulation) 
coders. In the first method, the ADM speech quality is shown to be 
improved by filtering the reconstructed (decoded) speech with a time- 
varying filter tailored to the natural speech bandwidth. In this case, 
adaptive bandpass filtering of the ADM output signal reduces coder 
noise by rejecting noise components at frequencies outside of the 
principal speech spectrum. Experimentally, we found that eliminating 
the upper 2 percent of the spectrum energy gave a reduction in average 
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bandwidth on the order of a factor of two relative to an initial 3-kHz 
bandwidth for typical speech samples. If the coder noise is white, there 
is an average noise power reduction by a factor proportional to the 
bandwidth reduction. Furthermore, the remaining portion of the noise 
power lies entirely within the band of the speech so that for reasonably 
good signal-to-noise ratios, some masking of the noise by the speech 
can be expected. 

The second case we ssilares is one of an adaptive sampling rate. In 
this case, the noise is again eliminated outside the principal speech 
bandwidth with a time-dependent low-pass filter. Then the average 
bit rate is reduced by a time-dependent decimation. 


1.1 Block diagram 


Figure 1 shows a block diagram of the system implemented in 
software for the tests to be presented. The system includes estimation 
of the short-time bandwidth (discussed in Section II), time-dependent 
filtering to this bandwidth, sampling rate conversion via decimation/ 
interpolation,’ and a 1-bit memory ADM coder with exponential step- . 
size adaption.” For. discursive purposes, we regard each of the two 
‘bandpass filters as a cascade of independently controlled low-pass and 
high-pass filters. The details on the implementation of this system are 
given in Appendix A. 


1.2 Test cases studied 
For clarity we define names for the following four cases studied: 


Normal ADM (ADM)—Both bandpass filters in Fig. 1 are fixed at the 
full voice-channel bandwidth, and the sampling rate is fixed. For 
example, 24 kbps ADM is implemented with a constant sampling rate 
and both filters are set to pass frequencies from 200 to 3200 Hz at all 
times (see Ref. 2). 

Post-filtered {ADM-PF)—Only the receiver reconstruction filter (at the 
far right in Fig. 1) varies to match the speech spectrum. The 
transmitter input filter is fixed at the channel bandwidth, and the 
sampling rate is fixed. | 

Pre- and Post-filtered (ADM-PPF)—Both the input and output filters 
are made to track the speech bandwidth, but -the sampling rate 
remains fixed as in ADM-PF. The addition of adaptive prefiltering 
allows the ADM coder to track the speech waveform with less error.° 

Adaptive Rate (ADM-AR)—The sampling rate varies at twice the upper 
cutoff frequency, and both the input and output filters track the 
speech bandwidth. 
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Fig. i—Block diagram of the variable-rate ADM coding system with adaptive bandpass filtering. The speech bandwidth is measured in real time 


to control the variable bandpass filters and ampling-rate conversion (decimation/interpolation). In the simplest case, only the right-most filter is 
varied along with the speech spectrum to provide coder noise suppression. In the case of adaptive sampling rate, both filters and the sampling rate 


are tailored to the speech. 


1.3 Results 


The main conclusions are: 

(t) Post-filtering gives a significant increase in quality. For example, 
16 kbps. ADM-PF gives a quality commensurate with 24 kbps apDM 
without post-filtering. The signal to noise ratio (s/n) is increased 
primarily in the low bandwidth segments such as back vowels, nasals, 
and voiced stops.* Almost all of the improvement arises from the 
adaptive low-pass component of the filtering. The adaptive high-pass 
filter contributes only slight noise suppression, and can introduce 
undesired side effects; for example, when there is a rapid transition 
from voiced to unvoiced, in which the speech band goes from low pass 
to high pass, an audible and objectionable change in the coder noise 
can occur even though there is an improved s/n due to the rejection of 
out-of-band low-frequency coder noise. Thus, adaptive low-pass filter- 
ing improves quality without serious side effects, while adaptive high- 
pass filtering only slightly improves s/n and causes significantly more 
audible noise modulation. For apM, these effects are most pronounced 
at 24 kbps and below. 

(it) For the case of ADM-AR (prefiltering, adaptive sampling rate, 
and post-filtering) we found that adapting the sampling rate causes 
low-pass time frames (such as voiced segments) to have a degraded 
s/n compared to the ADM s/n; however, for these frames, the bit rate 
is substantially reduced. Furthermore, while the in-band coder noise 
is increased, the out-of-band coder noise is eliminated. Consequently, 
the quality of ADM-AR is different from ADM but not easily judged to 
be worse. Informal listening tests indicated no reliable preference for 
one over the other for the few samples of speech tested (base bit rates 
of 16, 24, and 32 kbps). 

Summarizing positive practical results, our simulations indicate that 
time-dependent (adaptive) low-pass post-filtering yields a significant 
quality increase (for bit rates of 24 kbps and below) and adaptive low- 
pass pre- and post-filtering plus adaptive sampling rate yields reduced 
average bit rate with little change in quality. 

In Section II, we discuss how the time-dependent filter cutoff fre- 
quencies are measured from the short-time speech spectrum. Section 
III presents simulation results for the four cases defined above. 


ll. MEASUREMENT OF THE TIME-VARYING SPEECH BANDWIDTH 


Given the short-time spectrum of the speech at a given time, we 
wish to define the upper and lower cutoff frequencies of the spectrum 
in a way that minimizes bandwidth without introducing significant 
quality loss in the bandlimited speech. For this purpose, we define two 
constants 7’, and Ty which may be thought of as the fractional energy 
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of the speech bandwidth to be removed.’ T;, and Ty are taken to lie 
between 0 and 1. We call T, the lower cutoff threshold and Ty the 
upper cutoff threshold. If X(z, f) denotes the short-time spectrum of 
the speech at time ¢, then the time-varying high-pass and low-pass 
cutoff frequencies are found by solving 


aes ii 2 
Teas I. X(t, AI? af 


1 f(t) : 
TL -a5 | |X(t, f)|° df, (1) 


for fu(t) and fr(t), where E(t) is the total spectrum energy at time t 
given by 


E(t) -| [X(t f)[? df. (2) 
0 


Note that fu(t), the high-frequency (or low-pass) cutoff, and fr(t), the 
low-frequency cutoff, vary to maintain constant Tv, T. 

The discrete-time, discrete-frequency definitions that result from 
using the discrete Fourier transform (DFT) to generate short-time 
spectra are exactly analogous. When discussing sampled data, we write 
n in place of ¢, and the sampling rate will be denoted by f, = 1/T. 

The upper cutoff threshold Ty is the fixed fraction of the total 
energy that is rejected by the time-varying low-pass filter, and simi- 
larly, T; controls the time-varying high-pass filter. These constants 
are chosen in accord with desired coder quality. Ideally, the values of 
T., Tv might be optimized to trade off bandwidth for suppressed coder 
noise. In the spirit of Wiener filter theory, we might define the optimum 
thresholds as the values for which a decrease of either results in more 
added coder noise than added signal in the reconstructed speech, and 
where an increase of either value causes more distortion loss due to 
bandlimiting than quality gain from noise excision. However, we do 
not know how to define objective measures of subjective degradations 
due to changes in bandwidth and coder noise. In our tests, the thresh- 
olds T;, and Ty were set such that they did not appreciably degrade 
the speech quality in the absence of coder noise. That is, rather than 
attempt to define optimal thresholds for each bit rate, we wish merely 
to estimate the benefits of variable bandwidth when no perceptually 
significant distortion results from the bandlimiting alone. Accordingly, 
in all ADM coder simulations, where the initial speech bandwidth is 0.2 
to 3.2 kHz, the values T;, = 2 percent and Ty = 1 percent were used to 
specify the time-varying filters (and sampling rate when applicable). 

An example of the passband behavior for these threshold values is 
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given in Fig. 2. The phrase analyzed was from an adult male speaker. 
Note that within the telephone passband, the speech is basically either 
low pass or high pass at any given time. It is these temporal speech 
bandwidth variations that we exploit for noise and sampling rate 
reduction in ADM. 

Figure 3a gives a spectrogram of the same speech sample, and Fig. 
3b shows the spectrogram after filtering the speech to the bandlimits 
shown in Fig. 2. We see that the 1 percent energy upper cutoff limit 
tends to follow the third formant during voiced regions, and the 2 
percent lower cutoff limit has an almost unobservable effect. (The 2 
percent high pass has an audible effect on unvoiced phonemes, how- 
ever.) Figure 4a shows the output of a 16-kbps normal ADM coding 
system (time-invariant filters), and Fig. 4b shows the effect of post- 
filtering. The audible improvement due to post-filtering is much like 
Fig. 4b suggests, namely the out-of-band noise has been removed in 
Fig. 4b. This particular sample of post-filtered 16-kbps coded speech 
sounds about as good as when coded with normal 24-kbps ADM. 

Figure 5 gives a plot of the average bandlimits (averaged over the 
entire utterance) as a function of thresholds. When Ty = 
T, = 0, the passband is identical to the original speech passband; as 
the thresholds approach one half, the passband converges to zero. 
Comparing the two traces, we see that the speech is primarily low pass, 
which correlates with the fact that the utterance is predominantly 
voiced. Note that Fig. 5 implies an average bandwidth of only one half 
the maximum bandwidth using the values Ty = T, = 1 percent. In 
other words, it is possible to reduce the average sampling rate by a 
factor of two while sacrificing only 2 percent of the spectral energy. 

A few remarks are in order concerning practical issues associated 
with the measurement of the time-varying speech bandwidth. When 
tracking any spectrum over time, it is necessary to employ the proper 
balance of frequency resolution versus time resolution in the spectral 
analysis.* For speech, we wish to track bandwidth changes correspond- 
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Fig. 2—Spectral band edges vs time for a male speaker, obtained by rejecting the 
upper 1 percent and the lower 2 percent of the 200-3200 Hz spectrum energy. Band- 
edge values are computed every 12.5 ms. 
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Fig. 3—Spectrograms of 8 kHz sampled speech, before and after filtering to the 
measured speech bandwidth. (a) Original speech spectrum within the fixed channel 
bandwidth of 200-3200 Hz. (b) Same speech after filtering with a time-varying bandpass, 
which rejects the upper 1 percent and lower 2 percent of the spectrum energy in the 
band. This filtering is approximately transparent. 


ing to the articulation of phonemes. Tracking should be rapid so that 
there is little or no “smearing” of the estimated band-edges at the 
juncture of two dissimilar phonemes, and this implies using a small 
integration frame (window) in computing the short-time spectrum. 
However, when the spectrum is based on a frame length that is less 
than a pitch period, we obtain spectra that fluctuate excessively (at 
the pitch frequency) due to differing decay times of the vocal tract 
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Fig. 4—Spectrograms of the same speech as in Fig. 3 at the output of a 16 kbps apM 
system before and after time-varying filtering. (a) Fixed-bandwidth coder output. (b) 
Variable-bandwidth coder output. The time-varying filter is controlled by the same 
band edges as in Fig. 3b, i.e., the band edges are measured from the speech at the coder 
input. The effect of this filtering is to reduce significantly the ADM coder noise. 


impulse response components. For this reason, it is desirable to include 
at least 20 ms of speech in each spectrum computation, corresponding 
to the observation that the pitch of voiced speech rarely, if ever, falls 
below 50 Hz. As Fig. 6a shows, when the spectrum analysis integration 
time is less than a pitch period, the time-varying low pass cutoff fy(n) 
can oscillate quite significantly (e.g., +20 percent) at the pitch fre- 
quency. In Figs. 6b,c we show the effects of a seven-point median 
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smoother and a seven-point moving average smoother on the data of 
Fig. 6a. 

Another implementation issue is that of using the time-varying 
cutoff frequencies to control the output filters. Care must be taken to 
match the spectral analysis integration time, sampling interval for the 
cutoff frequencies, and filter impulse-response duration. These three 
times should be comparable in magnitude. In Ref. 6, it is shown that 
for the case of FFT-based (fast Fourier transform) analysis and filtering, 
the minimum sampling rate for the filter cutoff frequencies is deter- 
mined by the window used on the input to the FFT. For a length N FFr 
with a Hamming window, the band-edges may be sampled every N/4 
samples (i.e., successive FFTs used for calculating fu, f. may be offset 
in time by N/4 samples). It is shown in Ref. 7 that the resulting time- 
varying filter will have properly bandlimited coefficients regardless of 
the spectral modifications made on each FFT. 

Our formulation may be altered slightly to provide excellent per- 
formance during regions of silence. This is called the “idle channel” 
condition in the ADM literature. Inspection of (1) and (2) reveals that 
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Fig. 5—Time-averages of the upper and lower band edges vs spectrum energy- 
rejection threshold. The average is taken over the entire utterance of Fig. 2 for each 
threshold. At the far left of the figure, the energy-rejection thresholds are near zero so 
that the band edges lie at the outer extremes of the true speech band. At the far right, 
the two thresholds are equal to 0.5 corresponding to rejection of the upper and lower 50 
percent of the speech spectrum; consequently, the band edges meet at the median 
frequency. 
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Fig. 6—Ilustration of the effects of using an under-sized Fourier transform, and the 
possibility of compensation via filtering of the band-edge waveforms. (a) Upper 1 percent 
band edge for the same speech sample as in Fig. 2 using only a 10-ms time frame for 
spectrum computation. This causes oscillation of the measured band edge at the pitch 
frequency. (b) The same band-edge function of (a) after filtering with an order 7 median 
smoother. (c) The same curve of (a) after filtering with an order 7 moving average. 
Linear smoothing results in band-edge time behavior similar to that obtained when 
using a larger Fourier transform; however, the frequency resolution of the band-edge 
values is still sacrificed. 


the cutoff frequencies are indeterminate in this situation [|X (w, t)| = 
0]. If there is any amount of white noise present in the input signal, 
then the time-varying bandwidth, fu — fr, will open to full bandwidth 
as if the speech itself were spectrally flat. This undesirable behavior 
may be suppressed by various ad hoc schemes. In the simulations, we 
added a small positive value to F(t). That is, (2) is replaced by 


Bw =| X(t, f)|? df + omin, (3) 


0 


where oin may be thought of as noise energy, or as a lower bound on 
the acceptable speech level. As the speech energy falls to zero, the 
band edges cross, corresponding to disjoint low-pass and high-pass 
filters, and we must therefore define all negative values of fy — fi to be 
zero bandwidth. Also, there exists the possibility that no solution to 
(1) exists, for 02, > 0, in which case the bandwidth is again set to zero. 
Thus, when the channel is idle, there will be zero bandwidth and 
subsequently no output signal. This fact can be used to advantage 
when optimizing the step-size adaption algorithm in the ADM coder.® 
Our simulations assume that the band edges fy and f, are transmitted 
as side information in the variable bandwidth coding system. However, 
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it is worth considering filter adaption based only on the received 
speech data. Increasing the receiver energy rejection thresholds Ty 
and Ty, relative to the transmitter thresholds will contract the (esti- 
mated) band edges so as to compensate for the artificial band expan- 
sion that occurs because of coder noise in the received spectrum. 

If the bandlimits are transmitted as side information, then the 
increase in data rate is relatively small. As a practical example, if the 
FFT length is 512, a Hamming window is used, and the speech sampling 
rate is 8 kHz, then we have a pair of band edge values every 16 ms. 
Furthermore, the band edge values are quite smoothly behaved, and 
can be coded more efficiently. It appears that the band-edge waveform 
signals f(t) and fy(t) have a bandwidth on the order of 30 Hz for 
speech.° 


Ill. RESULTS OF ADM SIMULATIONS 


In this section, we present s/n evaluations of the four ADM coder 
configurations (described in Section I) ADM, ADM-PF, ADM-PPF, and 
ADM-AR. The comparisons are made using the s/n and segmental s/n® 
measures defined in Appendix B. Detailed parameter information may 
be found in Appendix C. 

The degree to which quality is enhanced by adaptive post-filtering 
depends on the character of the coder noise. If the coder noise is 
known to be stationary additive white noise, uncorrelated with the 
speech, then the gain in s/n may be predicted in advance from fy and 
ft. Given that bandlimiting the speech causes no distortion, the s/n of 
each segment will increase by 


(maximum bandwidth) 
(short-time bandwidth) 


= 10 log = 7 (4) 


where fmax is the full channel bandwidth. For a “typical” frame ( fu = 
2 kHz, fr = 400 Hz, C = 3 kHz), this is about 2 dB. The gain in quality 
at lower bit rates is dramatically greater than indicated by the s/n. 
This is perhaps due to the high perceptual significance of out-of-band 
coder noise and/or auditory masking of in-band noise by the speech. 

To anticipate the improvement of ADM due to post-filtering, we need 
to know the spectral distribution of ADM coder noise. While some 
theoretical work along these lines has been done,’ it is difficult analyt- 
ically to derive general estimates of the short-time noise power spectral 
density. Some intuition may be obtained, however, from simulations 
on isolated, quasi-stationary speech segments. 

Figure 7a shows a tenth-order Lpc spectral envelope” for the front 


s/n increase (dB) = 10 log 
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Fig. 7—Spectral envelopes of signal and noise for two representative sounds. All 
spectral envelopes were calculated using a tenth-order linear prediction on 1024 samples 
of data. The noise was isolated from the speech by subtracting the noiseless (precoder) 
speech from the apM coder output (which used time-invariant filtering). (a) Spectral 
envelope for the vowel “a” from “grab” superimposed with the spectral envelope of its 


669? 


associated ADM coder noise. (b) Spectral envelopes for the “‘s” sibilant in “sugar” and its 
corresponding coder noise. 


vowel “a” (as in “grab”) superimposed with the tenth-order spectral 
envelope of the normal 24-kbps apm coder noise generated by this 
vowel. (Appendix C gives detailed analysis parameters.) The measured 
s/n is 15 dB, and the noise spectrum within the passband is relatively 
flat. It should be noted that the slight ripple in the spectral envelope 
of the error signal depends on the order of the linear predictor. 
Figure 7b shows the same comparison of signal and noise spectral 
envelopes for the “sh” sound in “sugar.” Note that in this case, the 
noise is fairly flat out to 2.4 kHz after which it begins to follow the 
speech spectrum. The noise has a significant peak near 2.6 kHz 
indicating that these spectral components could not be properly 
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tracked. In this example, it was evident from the time-domain wave- 
form that the coder was tracking dominant high-frequency components 
with a large positive error in the amplitude difference estimate (adap- 
tive step-size inside the ADM coder).” The measured s/n for this sibilant 
is only 1.6 dB, and the primary character of the noise is that of rough 
loud “static.” 

Generalizing from Fig. 7, we might expect low-pass signals to gen- 
erate coder noise that may be approximately modeled as white, and 
high-pass speech segments to correspond to relatively strong correlated 
noise. Such heuristics, while over-simplified, serve to point out the 
more generally observed differences in ADM noise characteristics for 
voiced vs unvoiced speech. Awareness of these two contrasting cases 
aids in the interpretation of the segmental s/n in which the s/n for 
individual phonemes is evident. 

Figure 8b gives a plot of the segmental s/n (defined in Appendix B) 
versus time for the three cases ADM, ADM-PF, and ADM-AR. The bit 
rates of ADM and ADM-PF are 32 kbps. ADM-AR has 32 kbps as its 
maximum instantaneous bit-rate while the average rate for this partic- 
ular phrase is 23 kbps. Figure 8a shows the segmental input rms level 
from which the various phonemes may be located. The segment size is 
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Fig. 8—Time behavior of ADM noise for three cases as defined by the s/n of each 32- 
ms time frame (segmental s/n). (a) Segmental rms amplitude of speech utterance vs 
time indicating phoneme locations. (b) Segmental s/n vs time for normal ADM and post- 
filtered ADM (ADM-PF) at a bit rate of 32 kbps, and adaptive-rate ADM (ADM-AR) having 
a peak bit rate of 32 kbps. 
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32 ms in both Figs. 8a and 8b. We may observe several features in the 
behavior of the segmental s/n due to post-filtering and adaptive rate: 

(t) There is little difference among the three cases for front vowels 
such as “a” in “grab” and “e” at the beginning of “every.” From Fig. 
2, we see that during these segments, the speech bandwidth is wide 
and almost fully occupies the channel bandwidth. Consequently, the 
adaptive low pass is almost the same as the fixed low pass, and ADM- 
AR is running at maximum sampling rate during the greater portion of 
these vowels. 

(ti) When the low-pass cutoff fu(7) is small, ADM-PF realizes large 
quality gains due to rejection of much out-of-band coder noise. In 
contrast, ADM-AR exchanges these gains in return for reduced sampling 
rate. Examples of this may be seen at the phonemes corresponding to 
“hb,” “v,” “d,” and “u.” 

(tit) When the high-pass cutoff f,(7) is large [at which time fy(n) is 
maximum], ADM-AR reduces to the case ADM-PPF, and its performance 
is close to that of ADM-PF. Both exhibit higher segmental s/n than 
normal ADM due to elimination of low-frequency noise. This condition 
may be observed at the two unvoiced regions “sh” and “s.” 

We now turn to plots of segmental s/n averaged over the entire 
utterance, and we denote the average segmental s/n by S/Nseg. Figure 
9 shows S/MNseg VS bit rate for all four test cases. The post-filtered case, 
ADM-PF, is 2.8 dB better than normal ADM on the average. Note that 
the prefiltering in ADM-PPF, which reduces ADM tracking error, adds 
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Fig. 9—Segmental s/n averaged over the entire utterance for four cases, plotted as a 
function of bit rate. For the case of adaptive rate, the average bit rate is used as abscissa. 
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Fig. 10—Normal s/n computed on entire utterance for fsa cases, plotted as a function 
of bit rate. For the case of adaptive rate, the average bit rate is used as abscissa. 


still another dB or so to the S/Nseg for ADM-PF, and the improvement 
is always within the short-time speech band. The adaptive rate coder, 
ADM-AR, exhibits s/nseg between normal and post-filtered ADM. Overall, 
the s/n,g Measure corresponds well with subjective quality ratings. 
From informal listening to the speech samples represented in Fig. 9, 
we feel that ADM-PPF is not noticeably better than ADM-PF; ADM-PF is 
somewhat “cleaner” than ADM-AR (comparing where the maximum 
ADM-AR rate equals the ADM-PF rate), and normal ADM is definitely 
inferior due to the audible high-frequency noise which is allowed to 
pass. 

Figure 10 gives s/n (as opposed to S/Nseg) for the same four cases (cf. 
Appendix B). All deviations from Fig. 9 are due to the fact that the 
S/Nseg Measure is an average of the s/n’s (in dB) obtained from disjoint 
256-point frames while the s/n measure treats the entire speech sample 
as one frame. Since the regions of quality gain in ADM-PF and ADM-PPF 
are of low relative energy, they contribute little to the s/n. ADM-AR 
appears in this figure to be significantly superior, but this is misleading; 
ADM-AR has high distortion in the relatively low-energy low-bandwidth 
regions due to the large reduction in sampling rate, and the s/n 
measure does not adequately penalize it. For example, the consonants 
“b” and “d” might be least distinguishable in the ADM-AR case, relative 
to the other three, even though it scores the highest s/n. Thus, the s/ 
n measure is overly insensitive to low-amplitude intelligibility loss, 
especially in the case of ADM-AR. 
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IV. CONCLUSIONS 


It has been shown that the quality of ADM coded speech can be 
significantly improved by employing a time-dependent low-pass filter 
matched to the short-time speech bandwidth. A time-varying high- 
pass cutoff may be added with little additional computational cost, but 
its contribution to quality is small and sometimes perceptually dis- 
tracting due to audible noise modulation at bit rates below 24 kbps. 
Transmission of the slowly varying cutoff frequencies adds only slightly 
to the transmission bit rate. Two uses of the adaptive low-pass cutoff 
were discussed. First, time-varying low-pass filtering of the ADM de- 
coded signal was found to add quality commensurate with a large 
increase in ADM bit rate (e.g., 24 kbps quality at 16 kbps). Secondly, 
time-varying low-pass filtering before and after the ADM coder, coupled 
with a time-varying sampling rate, gave nearly the same quality as 
normal ADM but with a large reduction in the average bit rate (e.g., 24 
kbps from 32 kbps). The gains cited are for continuous speech, and 
better relative performance is to be expected for speech containing 
regions of silence. The final conclusions concerning quality are based 
on casual listening tests and are only indirectly supported by the s/n 
measures employed. 


V. ACKNOWLEDGMENTS 


The authors wish to thank N. S. Jayant for sharing his expertise on 
ADM coding, and for recurrent help throughout the project. We also 
thank J. L. Flanagan for numerous ideas and suggestions which were 
incorporated into this investigation. 


APPENDIX A 
Implementation of Variable Bandwidth ADM Simulation 


Referring again to Fig. 1, the software implementation is as follows. 
The input speech is sampled at 8 kHz, bandlimited to the typical 
channel bandwidth for telephone communication (200-3200 Hz), and 
is then resampled at 16, 24, 32, or 40 kHz. The data is partitioned into 
overlapping frames of 512 samples, a Hamming window is applied,” 
and the FFT of each frame is taken. The speech cutoff frequencies fu 
and f; are computed for each frame, as discussed in Section II. 

If prefiltering is included or if the sampling rate is to be lowered, the 
spectrum values outside the cutoff frequencies are tapered to zero 
using a precomputed filter band edge. The filter band edge is computed 
using a window design method based on a Kaiser window.” Next, an 
inverse FFT is taken on each frame, and the time-domain waveform is 
reconstructed by adding the frames back together, partially overlapped 
in time (overlap-add synthesis’). 
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The decimation stage is only active during sampling rate reduction, 
and it operates by selecting every mth sample, where m = [0.5 fs/fu] 
is the sampling rate reduction factor. [x] denotes the smallest integer 
=x, namely, m is the greatest integer such that m times the low-pass 
cutoff frequency for the current frame does not exceed the upper 
channel band edge. Note that the integer decimation method of varying 
the sampling rate does not take full advantage of the unused band- 
width; however, it has the advantage that it is quite simple to imple- 
ment. 

The coder is a one-bit ADM coder with exponential step-size adap- 
tation as described in Ref. 2. The coder output and the time-varying 
bandwidth information are assumed to be transmitted through a 
noiseless channel. 

The apm decoder is followed by a sample interpolator to restore the 
original sampling rate (when applicable), and the interpolator is fol- 
lowed by a time-varying filter. This filter is also implemented via 
short-time spectrum analysis, modification, and synthesis; it restricts 
the decoded speech spectrum to its original natural bandwidth, when 
post-filtering is employed, thus removing out-of-band coder noise. This 
filter is also part of the interpolation process as the interpolator merely 
inserts m = [0.5 f,/fu] zeros between each sample. 


APPENDIX B 
Signal-to-Noise Ratio Calculation 


Two types of signal-to-noise ratio are defined. The most common 
form is 


N-1 
XS (x(n) — px)? 

s/n 4 10 log as ; 
¥ (e(m) — pe)? 
m=0 


where x(m) is the signal with sample mean 
N-1 


1 
Lx = v2, x(m), 


e(m) is the noise with sample mean pe, and N is the total number of 
samples available for the s/n measurement. 

The definition of s/n diverges from subjective quality ratings for 
large N due to the fact that high-amplitude signal regions dominate 
the influence of low-amplitude signal regions during the s/n calcula- 
tion. This insensitivity may be partially circumvented by computing 
s/n values over segments of some reasonably small size M (e.g., 
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spanning 20 ms), and averaging the s/n (dB) values of the segments. 
Accordingly we define 


M-1 
Y [xe(m) — p(k)? 
segmental s/n(k) 4 10 log ——_—____—_ , 
Y Les) — wel)? 
m=0 
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where N is the number of segments of length M, x,(-) is the kth 
segment of the signal, and p(k) is the sample mean of the kth 
segment.’ In the case of S/Nseg, the measure is vulnerable to domination 
by segments having insignificant signal energy (i.e., the s/n can ap- 
proach —o in a time frame where the signal is silent and where there 
is any amount of noise). Consequently, if the total energy (sum of 
samples squared) in a given segment is below a prescribed energy 
threshold, the segment is eliminated from the computation of S/Ngeg. 
(This feature was not needed for the continuous speech samples used 
in the ADM simulations.) 

In all ADM tests, the noise e(m) is calculated as the point-wise 
difference between the noisy coded signal and a signal which was 
generated in precisely the same way but bypassing the ADM coder. In 
this way, all side effects of bandlimiting, processing delay, etc., are 
eliminated from the calculated error. Measurement of s/n in an ADM 
coding system is facilitated by the fact that it is a waveform coder (as 
opposed to source coder), and thus does not have the inherent delay, 
phase-dispersion, or level-offset characteristics that commonly impede 
the objective measurement of subjective signal quality. 


APPENDIX C 
System Parameters Used in Generating s/n Curves 


Coder input: Phrase = “Grab every dish of sugar” from an adult 
male speaker, sampled at 8 kHz, and bandlimited to 200-3200 Hz with 
a 256-point FIR bandpass. 

Time-varying filters: In all runs, the filters were implemented via 
modified FFTs of length N = 512. To prevent time-aliasing, the number 
of data points N, brought into the FFT input buffer plus the length N;, 
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of the Kaiser window (used as the basis of the time-varying filter) 
cannot exceed N. Furthermore, short-time spectral modification theory 
requires that the step-size through the data (time offset between 
successive FFTS) not exceed N,/4 for the case of a Hamming window 
on the FFT input.® The table below gives the employed data frame size 
N,, and time-varying filter length N, as a function of sampling rate f, 
for all ADM simulations. 


| fee) | Ne | 
| 16s] 804 208 
456 


112 








The filter controls fu(n) and fr(n) are each eight-bit values at a 
sampling rate of 4f,/ Nx. 

ADM coder: The step-size multipliers were experimentally found 
to give good results with P = 1.2, Q = 0.9.” These values did better 
than P = 1/Q = 1.5, P = 1/Q = 1.2, and a few other trial settings, in 
terms of the s/n and S/Nseg measures. 

LPC spectral envelopes: The short data segments and “s” 
were each processed with N = 512, f, = 24 kHz, N. = 456, Nz = 56, 
fixed filters, and nonadaptive sampling rate. The tenth-order LPc 
spectral envelopes were calculated using 1024 data samples. 


a” 
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One of the major drawbacks of the standard pattern-recognition 
approach to isolated word recognition is that poor performance is 
generally achieved for word vocabularies with acoustically similar 

_words. This poor performance is related to the pattern similarity 
(distance) algorithms that are generally used in which a global 
distance between the test pattern and each reference pattern is 
computed. Since acoustically similar words are, by definition, glob- 
ally similar, it is difficult to reliably discriminate such words, and a 
high error rate is obtained. By modifying the pattern-similarity 
algorithm so that the recognition decision is made in two passes, we 
can achieve improvements in discriminability among similar words. 
In particular, on the first pass the recognizer provides a set of global 
distance scores which are used to decide a class (or a set of possible 
classes) in which the spoken word is estimated to belong. On the 
second pass we use a locally weighted distance to provide optimal 
separation among words in the chosen class (or classes), and make 
the recognition decision on the basis of these local distance scores. 
For a highly complex vocabulary (letters of the alphabet, digits, and 
three command words), we obtain recognition improvements of from 
3 to 7 percent using the two-pass recognition strategy. 


1. INTRODUCTION 


As illustrated in Fig. 1, the “standard” pattern recognition approach 
to isolated word recognition is a three-step method consisting of 
feature measurement, pattern similarity determination, and a decision 
rule for choosing recognition candidates. This pattern recognition 
model has been applied to a wide variety of word recognition systems 
with great success.’* However, the simple, straightforward approach 
to word recognition, shown in Fig. 1, runs into difficulties for complex 
vocabularies, i.e., vocabularies with phonetically similar words. For 
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Fig. 1—Block diagram of standard approach to isolated word recognition. 


example, recognition of the vocabulary consisting of the letters of the 
alphabet would have problems with letters in the sets 


$1 = {A, J; K}, 

$2 = {B, C, D, E, G, P, V, T, Z}, 
os = {Q, U}, 

oa = {I, Y}, 

os = {L, M, N}, 

ge = {F, S, X}. 


Similarly, recognition of the computer terms of Gold® might lead to 
confusions among the set containing four, store, and core. In the above 
cases the problems are due to the inherent acoustic similarity (overlap) 
between sets of words in the vocabulary. It should be clear that this 
type of problem is essentially unrelated to vocabulary size (except 
when we approach very large vocabularies), since a large vocabulary 
may contain no similar words (e.g., the Japanese cities list of Itakura’), 
and a small vocabulary may contain many similar words (e.g., the 
letters of the alphabet). 

The purpose of this paper is to propose, discuss, and evaluate a 
modified approach to isolated word recognition in which a two-pass 
method is used. The output of the first recognition pass is an ordered 
set of word classes in which the unknown spoken word is estimated to 
have occurred, and the output of the second pass is an ordered list of 
word candidates within each class obtained from the first pass. The 
computation for the first pass is similar in nature but often reduced in 
magnitude from that required for the standard one-pass word recog- 
- nizer. The computation of the second pass consists of using an “opti- 
mally” determined word discriminator to separate words within the 
equivalence class. In Section H, we present the two-pass recognizer, 
and discuss its philosophy and method of implementation. In Section 
III, we give an evaluation of the effectiveness of the two-pass approach 
for a vocabulary consisting of the 26 letters of the alphabet, the 10 
digits, and the command words STOP, ERROR, and REPEAT. Finally, in 
Section IV, we summarize the results and show how they are applicable 
to practical speech recognition systems. q 
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ll. THE TWO-PASS RECOGNIZER 


Assume the word vocabulary consists of V words. The ith word, vj, 
is represented by the word template Ri, i = 1, 2, --- , V, where each R; 
is a multidimensional feature vector. Similarly, we denote the test 
pattern as T (corresponding to the spoken word q in the vocabulary), 
where T is again a multidimensional feature vector. For simplicity we 
assume that the pattern similarity and distance computation is carried 
out using the “normalize and warp” procedure described by Myers et 
al.,!° and illustrated in Fig. 2. A “standard” word duration of N frames 
is adopted, and each reference pattern is linearly warped to this 
duration. We call the warped reference patterns R;. Similarly, the test 
pattern is linearly warped to a duration of N frames, yielding the new 
pattern T. A dynamic time-warping alignment algorithm then com- 
putes the “standard” distance 


Soi | re 2 
D(T, Ri) = NW x d(T(k), Ri(w(k))) , (1) 


where d(T(z), R;(l)) is the local distance between frame k of the test 
pattern, and frame J of the ith reference pattern, and w(f) is the time- 
alignment mapping between frame k of the test pattern, and frame 
w(k) of the ith reference pattern. The total distance D of eq. (1) is 
only a function of 7. 

We define the local distance of the kth frame of the test pattern to 
the w(k)th frame of the ith reference pattern as d;(k), where 


d;(k) = d(T(k), Ri(w(k))), (2) 
so D(T, R;) of eq. (1) can be written as 
eS oo 1 
D(T, R;) = WV 2 d;(k). (3) 
N 






R,(n),n = 1,2,..., LINEAR 
TIME 


WARPING 








LINEAR 
TIME 
WARPING 









Tim), = 1,2,..,N 





Fig. 2—Block diagram of the normalize-and-warp procedure for equalizing the 
lengths of words. 


ISOLATED WORD RECOGNITION 741 


If R; corresponds to the correct reference for the spoken word T (i.e., 
i = q), then we would theoretically expect the local distance d,(k) to 
be independent of k, with d assuming values from a x” distribution 
with p (eight for the system we are using) degrees of freedom for the 
case where the speech features are those of an LPC model and the log 
likelihood distance measure is used for the local distance.” Thus, if 
we plotted d,(k) versus k, we would expect it to vary around some 
expected value d where 


d = E[d,(k)] = E[xp]. (4) 


An example of a typical curve of d,(k) versus k is given in Fig. 3a. 

If we now examine the typical behavior of the curve of d;(k) versus 
k when i ¥ q, we see that one of two types of behavior generally 
occurs. When word gq is acoustically very different from word 1, then 
di(k) is generally large [compared to d of eq. (4)] for all values of k, 
and the overall distance score D of eq. (3) is large. This case is 
illustrated in Fig. 3b. However, when we have acoustically similar 





k 


Fig. 3—Curves of d;(k) versus & for three cases. 
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words, then generally d;(k) will be approximately equal to d, (zk) for all 
values of k& in acoustically identical regions, and will be larger than 
d,(k) only in acoustically dissimilar regions. An example in which the 
dissimilar region occurs at the beginning of the word (the first N 
frames) is shown in Fig. 3c. 

The key point to be noted from the above discussion is that when 
the vocabulary contains words that are acoustically similar, and one of 
these similar words is spoken (i.e., it is the test utterance), then the 
total distance scores for these similar words consists of a random 
component [because of the variations of d(k) in the similar regions] 
and a deterministic difference (because of the differences in the dissim- 
ilar regions). In cases when the size of the dissimilar region is small 
(i.e., N < N in Fig. 3c), then the random component of the distance 
score can (and often does) outweigh the true difference component, 
causing a potential recognition error. For highly complex vocabularies 
(e.g., the letters of the alphabet), this situation occurs frequently. 

One possible solution to the above problem would be to modify the 
overall distance computation so that more weight is given to some 
regions of the pattern than others. For example, we could consider a 
weighted overall distance of the form 

N 

X WRk)d(T (Rk), Ri(w(k))) 

D(T, R;) = — (5) 

y Wie) 
k=1 
where W(k) is an arbitrary frame weighting function, and the denom- 
inator of eq. (5) is used for distance normalization. The problem with 
eq. (5) is that a “good” weighting function is difficult to define since 
the “optimal” set of weights is clearly a function of the “actually” 
spoken word (q) and the reference pattern being used (z). Furthermore, 
any weighting that would help discriminate between acoustically sim- 
ilar words, would tend to hurt the discrimination between acoustically 
different words. 

The above discussion suggests that a reasonable approach would be 
a two-pass recognition strategy in which the first pass would decide on 
an ordering of word “equivalence” classes (in which sets of acoustically 
similar words occurred), and the second pass would order the individ- 
ual words within each equivalence class. For the first-pass recognition 
an unweighted (normal) distance would be used, and for the second 
pass a weighted distance would be used. In order to implement such a 
two-pass recognizer, a number of important questions must be an- 
swered, including: 

(i) How do we “automatically” choose the word equivalence classes 
for each new vocabulary? 
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(ti) How do we determine class distance scores for the first recog- 
nition pass? 

(tit) How do we determine weighting functions for the second 
recognition pass? 

(iv) How do we generate weighted distance scores for the second 
recognition pass? 

(v) How do we combine results from both recognition passes to give 
a final, overall set of distance scores and word orderings? 

Some possible answers to each of these questions are given in the 
following sections. 


2.1 Generation of word equivalence classes 


Given the V vocabulary words v, v2, «++ , vy, we would like to find 
a procedure for mapping words into acoustic equivalence classes ¢,, 
j = 1, 2, ---, J, where J = V. There are at least two reasonable 
approaches for solving this problem; one is a theoretical approach, the 
other an experimental one. 

For the theoretical approach we can generate a “word-by-word” 
distance matrix D,,, on the basis of the phonetic transcriptions of the 
vocabulary entries. In order to do this we need to define a “phoneme” 
distance matrix, d,, a distance cost for inserting a phoneme, d;, and a 
distance cost for deleting a phoneme, dp. The phoneme distance matrix 
could be a count of the number of distinctive features that have to be 
changed to convert from one phoneme to another.” A total word-by- 
word distance is then defined by a dynamic time-warp match between 
the words, with a vertical step representing an insertion, and a hori- 
zontal step representing deletion. Figure 4a illustrates this procedure 
for the words eight and J, and Figure 4b for the words one and nine. 
For the words eight and J, the optimum path is an insertion (of -/), 
match between e’ and e’, and a deletion of ¢, giving a distance 


dr+ dp(e’, e’) + dp 


d(e't, Je’) = 3 , 


(6a) 


whereas for one and nine, the optimum path is a straight line giving 
dp(w, n) + dp(a, a’) + dp(n, n) 
3 : 


It should be clear that once d, (pi, pe), dr, and dp are defined, the word- 
by-word distance scores can be generated. 

A second approach to obtaining word-by-word distance scores is to 
use real tokens of the vocabulary words and do the actual dynamic 
time warping of the feature sets and obtain actual word distances. If 
several tokens have been recorded, averaging of distances increases 
the reliability of the final results. 


d(wan,na'n) = 


(6b) 
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Ap (el, e!) 


I, d 
. D dy + dp (el, ef) +dp 

dlett Je!) = Se ee 

3 
J dy 
el t 

(a) 

n dp{n, n) 
dp (W,n)+dp(a, al)ady(n.n) 

d(Wan, naln) = PO 

al dp (a, al) 3 





(b) 
Fig. 4—Examples illustrating “word” alignment based on dynamic time warping. 


From the word-by-word distance matrices, word equivalence classes 
may be obtained using the clustering procedures of Levinson et al.," 
in which the vocabulary words are grouped into clusters (equivalence 
sets) based entirely on pairwise distance scores. 

As an example of the use of the above techniques, consider the 39- 
word vocabulary consisting of the 26 letters of the alphabet, the 10 
digits, and the 3 command words STOP, ERROR, and REPEAT. These 39 
words become clustered into the sets 


Tokens 
oi = {B,C,D, BE, G, P, T, V, Z, 3, REPEAT}, 11 
go. = {A,J, K, 8, H}, 5 
63 = {F,S, X, 6}, 4 
go = {I, Bess 4}, 4 
5 = {Q, U, 2}, 3 
do = {L,M, N}, 3 
7 a {O}, 1 
os = {R}, 1 
oo — {W}, 1 
dio = {STOP}, 1 
oi. = {ERROR}, 1 
iz = {0}, 1 
gis = {1}, 1 
14 = {7}, 1 
dis = {9}. 1 
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We discuss this vocabulary and the resulting equivalence sets a great 
deal more in Section III. 


2.2 Determination of class-distance scores 


Once all the vocabulary words have been assigned to one of the J 
classes, the first recognition pass estimates an ordering of the word 
classes in terms of class-distance scores. The class-distance scores can 
be determined in one of two ways. First they can be computed as the 
minimum of the word-distance scores, for all words in the class, i.e., 

d(¢;)= min D(T,R;), j=1,2,-+- oa. (7) 
U; E $; 

This computation is similar to the one used by Aldefeld et al.” for 
directory listing retrieval. 

An alternative method of obtaining class-distance scores would be 
to obtain “class-reference” templates (as well as word-reference tem- 
plates) and to measure distance directly from the class-reference 
templates. Clearly with multiple templates per class, the K-nearest 
neighbor (KNN) rule can be used as effectively for class templates as 
for word templates. 

The reason for considering class-reference templates for obtaining 
the class-distance scores is that the number of word classes is clearly 
smaller than the number of words. Hence, the number of distance 
calculations required to establish class-distance scores is generally 
much lower for class templates than for word templates. For example, 
for the 39-word vocabulary discussed previously, there are 15 word 
classes. Hence there is almost a 3 to 1 reduction from words to word 
classes. However, it should be clear that the danger in using class 
templates is that errors in determining class distances can be made 
from the reduced number of templates. This point will be discussed 
later in this paper. 


2.3 Choice of weighting functions for the second pass of recognition 


The output of the first recognition pass is an ordered set of word 
class-distance scores. For the second recognition pass, all words within 
the top class (or classes) are compared to the unknown test-word 
pattern (T) using a weighted distance of the type discussed in eq. (5), 
and an ordering of words within the class is made. If several classes 
have similar class distance scores, the words within each of these 
classes are ordered in the same manner. 

_ The key question that remains is how do we choose the weighting 
function, W(k), of eq. (5) in an optimal or reasonable manner. The 
reader should recall, at this point, that the optimal weighting function, 
W(k), is assumed to be a function of the pair of indices z (the reference 
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pidf=;) pidfi#)) 





Fig. 5—Simple Gaussian model for frame distance distributions. 


word) and 7 (the proposed test word). Hence if there are L words in an 
equivalence class, then there are L(LZ — 1) sets of weighting functions 
{the cases 1 = j have W(k) = 1]. 

We have investigated two ways of determining W(k) for the second 
recognition pass. Optimality theory says that to maximize the weighted 
distance of eq. (5),'° the value of W(k) should be 


W(k) =1 k= Ro, (8a) 
=0 all other R, (8b) 


where ko is the index where the distance between R; and T is, on 
average, the maximum. In this manner, the algorithm places all its 
reliance on the single frame where one would expect the maximum 
difference between reference and test patterns to occur. In practice, 
this weighting does not work since the variability in location of the 
frame k = ky of eq. (7) is large. Hence, on several trials the distances, 
using the weighting of eq. (7), can vary considerably. 

A more effective manner of determining a good (but not optimal) 
set of weights is as follows. Consider the model for the distribution of 
distances for a single frame as shown in Fig. 5. The curve on the left 
in Fig. 5 is the assumed distribution of distances in the case when 
i=] (ie., the reference and test patterns are from the same word). In 
this case, we expect a x” distribution with p (order of the Lpc model) 
degrees of freedom for the frame distance. For convenience, we model 
this distribution as a Gaussian distribution with mean m, and standard 
deviation 0.” 

For the case when 1 +7 (i.e., the reference and test patterns are from 
different words), we assume the frame distance has a Gaussian distri- 


* This assumption is reasonable since the word distance, which is a sum of frame 
distances, has a Gaussian distribution (by the central limit theorem), and the actual 
probability of word error is directly related to the word distance. 
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bution (as shown to the right in Fig. 5) with mean mz and standard 
deviation oz. 

We now make a simple recognition model that says the probability 
of recognition error for the word is proportional to the probability of 
error for single frames (since the word distance is the sum of frame 
distances). Then, based on the model of Fig. 5 with assumed Gaussian 
statistics, the probability of correct classification (i.e., finding a smaller 
frame distance for the spoken word, than for any other word) for a 
single frame is 


P(C) = | P[p(d/inj) = A]- Pl p(d/inj) > A] dd, (9) 


where P[x] is the probability of the event x occurring. Equation (9) 
says that the probability of correct frame classification is the integral 
of the probability that for the correct word (i = 7) we get a frame 
distance A, and for the closest incorrect word (i ~ 7) we get a frame 
distance greater than A. Thus the probability of a frame error is 


P(E) =1-P(C), (10) 


which becomes 
P(E) = | N[A — mm, a | N[n — me, 62] dyn dX, (11) 
—oo A 


which can be put into the form 


(m2—m,) /(oj+03)'? 2 : 
exp (—x°/2) me —- mM, 
P(E) -| —_——_ dx = Erf |= }.. (12) 
sae V27 Vo? + 03 


The form of eq. (12) can be verified for the simple cases mz = m,, 
where P(E) = 0.5, and mz >> m, where P(E) —> 0. 

The above discussion suggests that a reasonable choice for frame 
weighting would be 


| (dii(k)) — (di(R)) | 


(04,,(2) + 04,(%))” 


Ww'(k) = (13) 


where dj;(k) is the local distance between repetitions of word 7 for 
frame k, and dij (k) is the local distance between spoken words / and i 
for frame k, and where the expectations are performed statistically 
over a large number of occurrences of the words vu; and U;. 

By way of example, Fig. 6 shows examples of plots of (dj:(R)) versus 
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DISTANCE 





DISTANCE 





-- A, J, K, &8vs A 


DISTANCE 





FRAME NUMBER (k) 


Fig. 6—Examples of frame-by-frame distances for words within word equivalence 
classes. 


k and W(k) versus k for some typical cases.* Figure 6 shows a series 
of plots for the following cases: 


(i) (Fig. 6a) Curves of (dji(k)) and od) for the case where word L 
was the letter J, and word j was the letter Y. We can see that (dii(R)) 
(the solid curves) is approximately constant whereas (dji(R)) differs 
from (dji(k)) only at the beginning of the word (i.e., the first eight 


* The data of Fig. 6 were obtained from about 10,000 comparisons for each word, i.e., 
a large data base was used. 


ISOLATED WORD RECOGNITION 749 


frames). We also see that the curves of 04,2) (the dashed curves) are 
comparable for the cases j = 7 and for 7 ¥ 1, with only small differences 
occurring in the first eight frames. 

(ii) (Fig. 6b) Curves of (dj;(k)) and O9,(h) for the case where word i 
was the letter A, and where j corresponded to the letters J and K for 
word 8. Similar behavior to that of Fig. 6a is seen, in that (dji(#)) is 
approximately constant, and (dj;(k)) is larger than (dj;(k)) at the 
beginning of the word, for words J and K, and at the end of the word, 
for word 8. For the word 8, the curve of o4,;z) is also fairly large at the 
end of the word, indicating the high degree of variability in the plosive 
release of the word 8. 

(tii) (Fig. 6c) The part shows the results of averaging the data of 
Fig. 6b over all 7 + 1 with in the class of word i, i.e., class-weighting 
templates. In this case the curve of (dj;(k)) shows flat behavior except 
at the beginning (due to J, K) and end (due to 8). If storage of word- 
weighting curves is burdensome, the use of class-weighting curves 
could be considered as a viable alternative. 

Figure 7 shows a set of two weighting curves W”'(k) for the words 
I and Y. Figure 7a shows the weighting curve for reference word J and 
test word Y, and Fig. 7b shows the weighting curve for reference word 
Y and test word J. Several interesting properties of the curves should 
be noted. First we see that W”'(k) generally consists of a large pulse 
(for these examples this occurs near k = 1) and a residual tail. The tail 
is a measure of the statistical noise level, i.e., the statistical difference 
between (dj:(k)) and (dj;(k)) in the region of acoustical similarity. 
Typically the peak amplitude in the tails is less than 10 percent of the 
peak amplitude in the main pulse. 

Another interesting property of the weighting curves is that there is 
no symmetry, in that 


W*/(k) ¥ W"(R). (14) 


An explanation of this behavior is given in Fig. 8, which shows two 
plots of dynamic time-warping paths for the words J and Y, where it 
is assumed that the word Y is simply the word J with a prefix phoneme 
/w/. Figure 8a shows that when J is warped to Y, there is a discrepancy 
region in which the /w/ is being warped to the initial region of the 
/a'/ and large distances result. The /a7/ is warped to itself (the 
“ideal” path) and no further distance is accumulated. Figure 8b shows 
that the discrepancy region is considerably smaller when Y is mapped 
to I. The resulting weighting curves agree in form with the results 
given in Fig. 7. 


2.4 Generation of distance scores for the second recognition pass 


We have now shown how to assign words to classes, how to get class 
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WEIGHT 


REFERENCE: Y 
TEST: J 


WEIGHT 





FRAME NUMBER (k) 


Fig. 7—Weighting curves for comparing the words /I/ and / Y/. 


distance scores for the first recognition pass, and how to assign weights 
for pairs of words within a word class. The next step in the procedure 
is the determination of the distance for the second recognition pass 
based on the pairwise weighted distance scores. 
To see how this is accomplished, we define a pairwise weighted 
distance Dj; as 
N 


» W""(k)di(k) 


k=1 


Di = : (15) 


N so 
YW? (Rk) 
k=1 

where 7 is the index of the reference pattern (i.e., one of the words in 
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Fig. 8—An example showing why the word weighting curves are not symmetrical. 


the equivalence class) and / is the (assumed 1) ind of He Test) natant 
(again one of the words in the equivalence. class). _- 

The quantity D;; of eq. ( (15)-is Comput ed for all i, 7 pairs (with i ¥/) 
in the wo d class with riinimum class distance, and a matrix of pairwise 
oo is obtained. The word distance, D,, can be obtained in one 
of two ways, namely: 

(tz) Averaging over the 7 index, giving 


D; = > Di (16a) 


JAt 


(ti) Finding the minimum over the 7 index, i.e., 


D; = min {D,,i}. (16b) 
J 


JA 


The advantage of averaging is that D; tends to be more reliable, since 
averaging is equivalent to adding weighted distances over a larger 
number of frames than would be used for a single comparison. The 
minimum computation is useful, especially when several of the D,; are 
about the same. We examine oe these scoring methods in Section 
Til. 

For the case of averaging pairwise distance scores [eq. (16a)], the 
computation can be carried out more efficiently as follows. By combin- 
ing eqs. (15) and (16a) we get 
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N 
> W*(k) di(k) 
k=1 


D; = » D;i = » a ee (17a) 
J J y W?(k) 
k=1 
N (We (k) di(k 
Ss dda RNa) (17b) 
J k=1 5 Wi (k) 
k=1 
N Wi(R 
=e) ENG (17c) 
mT NY WH) 
k=1 
N et ie 
= VF Wk) dik), (17d) 
k=1 
where 
ant wr 
Wieoe (18) 
’ SW (Rk) 
k=) 


Thus, for L words in the equivalence class, we can compute D; with N 
multiplications and additions [rather than the N(L — 1) computations 
of eq. (16a)], and only L vectors of N averaged weights [W'(k)] need 
be stored, rather than L(L — 1) vectors as implied by eq. (15). 

Another variation on the distance weighting that was studied here 
was the effect of applying a nonlinearity to the weighting function, 
W*", before computing D,,;. The nonlinearity was to replace W”(k) by 
W’*"(k), defined as 


bres ji if W""(k)/W, > T, 
WHR) = i (k) i ( )/ Wax (19) 
0 otherwise, 7 
where 
Wmax = max [ W*"(k)], (20) 


and T is a threshold which is specified in the algorithm. The nonline- 
arity of eq. (19) truncates (to 0) the weighting curve whenever its 
relative amplitude falls below the threshold. Figure 9 illustrates a 
typical curve W*"'(k) and its truncated version W”'(k). The new 
weighting function was then applied directly in eq. (15) in place of 
W’' (k). Clearly, when T = 0, W”'(k) and W*"(k) are identical. Again, 
when averaging is used, the computation of eq. (17) gives a reduced 
set of weights. 
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Fig. 9—An example of a weighting curve and its truncated version. 


2.5 Overall distance computation 


If we can make the assumption that the probability of a class error 
on the first recognition pass is significantly smaller than the probability 
of a word error on the first pass, then the final distance for each word 
of the minimum class is the distance obtained on the second recogni- 
tion pass. However there are applications in which it is desirable to 
have a distance score for every word in the vocabulary. Hence, in these 
cases, it is necessary to combine the ordering from the second pass, 
with the distances from the first pass. The basis for such a strategy is 
that distances on the first pass are statistically more reliable than 
distances on the second pass, whereas order statistics (within the class) 
are more reliable on the second pass than on the first pass. One very 
simple way of combining distances and word orders is to obtain second- 
pass ordering for every word in the vocabulary (i.e., apply the method 
of Section 2.4 to all word classes), and then reorder the word list using 
distances from the first pass, and ordering within the class from the 
second pass. 


2.6 An example of the use of the two-pass system 


To illustrate this entire procedure, Tables I to [II show an example 
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Table |—Recognition results for a simple example 
(first pass) 


Word 
Word Position Class 
Word Word Distance First Class Distance 
Index Class First Pass Pass Number First Pass 


0.47 1 0.47 
0.39 2 0.66 
0.51 3 0.37 
0.72 
0.42 
0.60 
0.67 
0.83 
0.37 
0.78 
0.66 
0.62 


— 


COON WHE 
_ 
NWR NOMDWOUON 


10 
12 


MSMbDnw WM rR tN wwe 
— 


of the recognition steps for a 12-word vocabulary with three word 
equivalence classes. Table I shows the results of the first recognition 
pass. The class distance scores are assigned as the minimum word 
distance for words within the class. The “best” class in the first pass 
is class 3 with a distance score of 0.37, with class 2 having a somewhat 
higher distance of 0.47. In the second recognition pass the words within 
the best class (or classes) are compared using the optimally determined 


Table II—Second recognition pass results for the example in Table | 





1 
i | 0.60 3 Class 1 
| 0.69 | 4 
| 0.60, 2 
Dj D;(avg) Order 


3 
l 2 Class 2 
4 
1 
Order 
1 
i 4 Class 3 
3 
2 
Order 
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weighting functions. The results for each of the three classes are shown 
in Table II. In practice, one would usually need to compute the D;,; 
scores only for the best one or two classes. However, for explanatory 
purposes, results are shown for all three classes. Also, as discussed 
above, in the case of distance averaging, the D;; scores need not be 
computed since the D; scores can be obtained directly via eqs. (17) and 
(18). Using the technique of averaging leads to the within-class dis- 
tances and orderings as shown in the table. Finally, Table III shows 
the results of reordering the words using the distances obtained from 
pass 1, and the within-class orderings obtained from pass 2. Thus word 
2 is the best recognition candidate (with a distance of 0.37), whereas 
word 9 was the best recognition candidate at the end of the first pass. 
Other, within-class reshufflings of word position occur as a result of 
the two recognition passes as shown in Table I. 


2.7 Summary of the two-pass recognizer 


Figure 10 shows a block diagram of the full two-pass isolated word 
recognition system. In the first pass a DTW distance is computed 
between the unknown test word and the reference templates for each 
word class. The outputs of the first pass are ordered sets of word 
distance scores and class distance scores. 

For the second pass a set of pairwise weighted distances is deter- 
mined for all words within each word class with suitably low scores on 
the first recognition pass. The final recognition output is a combination 
of distance scores from the first pass and word orderings from the 
second pass. In the next section we demonstrate how this procedure 
works in some practical recognition examples. 


Table II!—Overall word 
positions and distances for 
the example given in Tables 

| and Il 


Word Word Word 
Index Position Distance 


0.47 
0.37 
0.51 
0.78 
0.42 
0.62 
0.67 
0.72 
0.39 
0.83 
0.66 
0.60 


ro 
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Fig. 10—Block diagram of the overall two-pass recognizer. 


il. EVALUATION OF THE TWO-PASS RECOGNIZER 

To test the ideas behind the two-pass recognizer, we used a data 
base of existing recordings. The word vocabulary consisted of the 
V = 39 word vocabulary of the letters of the alphabet, the digits (0 to 
9), and the three command words STOP, ERROR, and REPEAT. The 
training data for obtaining word and class reference templates, and 
pairwise word weighting curves, consisted of one replication of each 
word by each of 100 talkers (50 men, 50 women).* The word reference 
templates (12 per word) were obtained from a clustering analysis of 
the training data.'** A set of “class” reference templates (12 per class) 
was obtained from a second clustering analysis in which the words 
within a class were combined prior to the clustering. The pairwise 
word weighting curves were obtained by cross-comparing all word 
tokens within a word class, averaging the time-aligned distance curves, 
and computing both the averages and standard deviations for each 
frame. 

To test the performance of the overall system, two test sets of data 
were used. These included: 

1. TsS1—10 talkers (not used in the training) spoke the vocabulary 
one time over a dialed-up telephone line. 

2. TS2—10 talkers (included in the training) spoke the vocabulary 
one time over a dialed-up telephone line. 

Two sets of performance statistics were measured. For the first 
recognition pass the ability of the recognizer to determine the correct 
word class was measured. For the second recognition pass the improve- 
ment in word recognition accuracy (over the standard one-pass ap- 
proach) was measured. The results obtained are presented in the next 
two sections. 


* All results presented here are for speaker independent systems. 
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Fig. 11—Plots of class accuracy as a function of the number of templates per word 
(Q), class position (C), and KNN rule (KNN) for a 15-class vocabulary. 


3.1 Class recognition accuracy for the first pass 


The ability of the recognizer to determine the “correct” word class 
of the spoken word was measured using both word templates (and 
obtaining class-distance scores from the word-distance scores as dis- 
cussed previously), and class templates (obtaining class-distance scores 
directly). The number of templates per word (or per class) varied from 
1 to 12 in the tests to see the effects of the number of reference 
templates on the class accuracy. The K-nearest neighbor (KNN) rule 
was used to measure class scores with values of KNN = 1 (minimum 
distance), KNN = 2 (average of two best scores), and KNN = Q (average 
of Q best scores), where Q was the total number of templates used per 
word (or per class). 

The results of the class recognition accuracy tests are given in Figs. 
11 and 12.* Figures 11 and 12 show plots of class error rate (based on 
the top C classes) as a function of the number of templates per word 
(Fig. 11) or templates per class (Fig. 12), for values of KNN = 1 and 2, 
and for C = 1 (top candidate), C = 2 (two best classes), and C = 3 
(three best classes). Figure 11 shows results when each class is repre- 
sented by word templates, and Fig. 12 shows results when each class 
is represented by class templates. 


* The reader should note the difference in vertical scales between Figs. 11 and 12. 
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Several interesting observations can be made from Figs. 11 and 12. 
These include: 

(t) The KNN = 1 rule performs consistently better than the KNN = 
2 rule for class discrimination, for all values of C and Q. This result is 
in contradiction with the results of Rabiner et al.° who found signifi- 
cantly better performance for KNN = 2 than for KNN = 1. The 
explanation of this behavior is that the KNN = 2 rule provides signifi- 
cantly improved, within-class discrimination (at the expense of slightly 
worse between class discrimination), and that when the only function 
is to determine the class, the KNN = 1 rule is superior. In fact when the 
KNN rule was used with a value of KNN = Q (i.e., averaging over all Q 
reference templates), the class accuracy on the first candidate de- 
creased by about 20 percent—a highly significant loss of accuracy. 
This result again demonstrates that the minimum distance rule 
(KNN = 1) is best for class discrimination. 

(ii) The use of word-reference templates provides significantly better 
performance than obtained from class-reference templates. For ex- 
ample, the class error rate for the top three classes (C = 3) with Q = 
4 templates per word is essentially 0; whereas the class error rate for 
the top three classes with four templates per class is about 4 percent. 
This result shows clearly the importance of representing each word in 


CLASS ERROR RATE IN PERCENT 





NUMBER OF TEMPLATES PER CLASS (Q) 


Fig. 12—Plots of class accuracy as a function of the number of templates per class 
(Q), class position (C), and KNN rule (KNN) for a 15-class vocabulary. 
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the equivalence class by an adequate number of word-reference tem- 
plates. 

(tii) With six templates per word, error rates of about 4 percent 
(C = 1), 1 percent (C = 2), and 0 percent (C = 3) are obtainable, 
indicating that the full contingent of 12 templates per word is unnec- 
essary for proper class determination. Using 6, rather than 12 templates 
per word reduces the computation in the first recognition pass by 50 
percent. If we always use two or more word classes, the required 
number of templates per word for the first pass can be reduced to four, 
with no serious loss in class accuracy. 

The results shown in Fig. 11 indicate that high accuracy can readily 
be achieved in determining the correct equivalence class for each word 
in a very complex vocabulary. Hence there would appear to be no 
problems in implementing the first pass of the recognition system. 


3.2 Within-class word discrimination for the second pass and overall 
performance scores 

The two-pass word recognizer was tested on the words of Ts1 and 
TS2. For each test set a total of 390 words were used (39 words X 10 
talkers). For TS1, the word recognition accuracy (for the best candi- 
date) on the first pass was 78 percent, and for Ts2 (with talkers from 
the training set) the word recognition accuracy on the first pass was 85 
percent. At the output of the second pass, the word recognition 
accuracy for the best candidate [using the averaging technique of eq. 
(16a) and assuming the correct word equivalence class was found] was 
84.6 percent for Tsl and 88.5 percent for TS2, representing potential 
improvements of 6.6 percent and 3.5 percent, respectively. The reason 
that a larger improvement in accuracy was obtained for Ts1 data than 
for TS2 data was that the accuracy on the first pass was lower for Ts1 
than for TS2 (where the talkers were in the training set) and hence 
there was more room for improvement within the word classes. 

Figures 13 and 14 show plots of the changes in accuracy that are 
obtained for Ts1 (Fig. 13) and Ts2 (Fig. 14) data when a threshold is 
imposed on the distance scores at the output of the first recognition 
pass. The threshold specifies that the second recognition pass is 
skipped if the distance of the second word candidate is more than the 
threshold greater than the distance of the first word candidate. Clearly 
this procedure is a strictly computational one, since low-distance scores 
for a single word on the first pass are highly reliable indicators that no 
second pass is necessary. The data plotted in Figs. 13 and 14 show the 
percentage of cases where the actual spoken word comes in a lower 
position on the second pass than in the first pass within the word class; 
it also shows the percentage of cases when the spoken word comes in 
a higher position on the second pass than the first pass, and the 
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Fig. 13—Percentage improvement, decrease, and the resulting difference in word 
position at the output of the second recognition pass for TSs1 data as a function of the 
distance threshold using the averaging method. 


difference (the improvement) between the two curves. All the results 
are plotted as a function of the distance threshold for performing the 
second-pass computation. It can be seen from these figures that the 
two-pass recognizer is not ideal, i.e., there is a significant fraction of 
words for which a worse position results at the output of the second 
pass. However, on balance, it is seen that a real improvement in 
recognition accuracy results, and it is this improvement that makes 
the procedure a viable one. 

A similar set of results obtained using the minimum computation of 
eq. (16b) on the second pass rather than the average computation of 
eq. (16a) are shown in Figs. 15 and 16 for Tsl and Ts2, respectively. 
These plots show the same information as those of Figs. 13 and 14 for 
the averaging procedure. A comparison of these results shows that the 
averaging computation performs as well as, or better than, the mini- 
mum computation for the whole range of distance thresholds, and for 
both data sets. These results indicate that the averaging method 
provides a small but important statistical stability to the computation. 


3.3 The effect of thresholding on the weighting curves 


We ran a series of tests with the data from Ts1 and Ts2 to investigate 
the effects of applying thresholds to the weighting curves as illustrated 
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Fig. 14—Percentage improvement, decrease, and the resulting difference in word 
position at the output of the second recognition pass for TS2 data as a function of the 
distance threshold using the averaging method. 


in Fig. 9. The results indicated that poorer performance always re- 
sulted when any significant part of the weighting curve was zeroed out. 
Thus the gain achieved by removing the “statistical” low-level parts of 
the weighting curve was canceled by the “deterministic” loss from the 
rest of the weighting curve. Hence the conclusion was to use the entire 
weighting curve as derived from the statistical model. 


3.4 Computation for the two-pass recognizer 


We have seen in Section 3.3 that word recognition accuracy improve- 
ments of from 3.5 to 6.6 percent result for the 39-word vocabulary 
using the two-pass recognizer. A key question that must be answered 
is what is the cost of the computation for the two-pass system. 

To answer this question we must examine the computation in each 
pass of the recognizer. In the first recognition pass, for a V-word 
vocabulary with Q templates per word, a total of QV ptw comparisons 
are made. For a value of N = 40, each DTW comparison requires about 
500 nine-point dot-product computations, so a total rate, Ri, of 


R, = Q-V-500-9 (21) 


multiplications and additions are required. 
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If we assume that the local distances d;;(k) associated with the 
optimum warping paths are saved for each reference template, then 
for each pairwise comparison of the second pass a total of N (typically 
40) multiplications and additions are required. For L words in the 
equivalence class, a total of 


R,=L-(L—-—1)-N (22) 


multiplications and additions are required for the second-pass com- 
putation for a single equivalence class. For the averaging procedure of 
eq. (17), Re is reduced to LN multiplications and additions. 

If we assume typical values of V = 39, @ = 12, L = 7, N = 40, we get 
R, = 2,106,000 and Rz = 1680, i.e., the computation of the second pass 
is insignificant compared to the first pass computation. Furthermore 
since we can use reduced values of Q for the first pass (i.e., Q = 6 or 
Q = 4) the overall computation can be significantly reduced from the 
standard isolated word recognizer, with the same improvement in 
accuracy! 


IV. DISCUSSION 


The results presented in the preceding section show that improved 
recognition accuracy can be obtained via a two-pass recognition algo- 
rithm. It was shown that the improvements were both global, i.e., in 
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Fig. 15—The same results as in Fig. 13 obtained using the minimum method. 
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Fig. 16—The same results as in Fig. 14 obtained using the minimum method. 


an absolute recognition sense, and local, i.e., within the classes of 
equivalent words. Although the proposed two-pass recognizer has a 
number of possible implementations, it was shown that the best choices 
were to use a reduced set of word templates on the first pass, and to 
use all word classes that had reasonably small distance scores on the 
second pass. 

One of the major issues that remains unresolved in the two-pass 
recognizer is the choice of weighting curve used in the second-pass 
distance computation. The assumed Gaussian model which led to the 
variance-weighted difference of means for the weights is, at best, an 
‘approximation to the actual situation. Experimentation with modified 
forms of the weighting curve of eq. (13) led to poorer recognition 
performance. Thus, because we lacked a viable alternative, the weight- 
ing curve of eq. (13) is the only one we investigated for use in the two- 
pass recognizer. 

An interesting question that arises as a result of this study is how 
could this two-pass recognizer aid in practical recognition tasks. As 
one would anticipate, the answer to this question is that it depends on 
the specific recognition task. For example, for the backtracking direc- 
tory listing retrieval system of Rosenberg and Schmidt,” the improve- 
ment in recognition accuracy could provide significant reductions in 
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search time. However, for the search procedure of Aldefeld et al.,’° the 
increased word accuracy would have no effect on the search time, but 
could increase the name accuracy, especially when similar names exist 
in the directory (e.g., T. Smith and P. Smith). For applications like the 
airlines reservation system of Levinson and Rosenberg,’ the increased 
word accuracy would reduce the load on the syntax analyzer; however, 
it needn’t necessarily increase the overall accuracy of the system. 

The above examples show that the two-pass recognition strategy 
can be useful for some applications, but one must examine carefully 
the specific task before claiming how useful it will potentially be. 


V. SUMMARY 


We have shown that a two-pass approach to isolated word recogni- 
tion is viable when the word vocabulary consists of sets of acoustically 
similar words. The first recognition pass attempts to determine accu- 
rately the class within which the spoken word occurs, and the second 
recognition pass attempts to order the words within the class, based 
on weighted distances of pairwise comparisons of all words within the 
class. We discussed several alternatives for implementing this two-pass 
recognizer, and we made a performance evaluation which showed that 
a reliable class decision could be made based on a reduced set of 
template scores, and an improved word decision could be made from 
weighted pairwise distance scores. 
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l. INTRODUCTION 


Public-key cryptosystems offer a degree of flexibility not available 
with conventional (private-key) systems.’” In particular, the key re- 
quired for decryption in a public-key system can be changed at will, 
even in the middle of a message. This makes the task of the eaves- 
dropper very difficult indeed. A frequently cited disadvantage of pub- 
lic-key systems is their relative slowness (typically a few kilobit/s) 
caused by the large amount of number-crunching they require.** This 
has led to the development of hybrid cryptosystems in which a key, 
exchanged via a slow public-key system, is subsequently used in a fast 
conventional system, such as the Data Encryption Standard (DEs).” In 
this paper we present a fast algorithm for executing the knapsack 
cipher (a public-key cryptosystem).° When implemented with TTL 
integrated circuitry, this algorithm should permit data rates in the 
neighborhood of 10 Mbit/s. This speed is sufficient to provide security 
for a wide range of voice, data, and narrowband video traffic without 
the need for a hybrid cryptosystem. 

Section II presents an elementary example of the knapsack cipher 
to show how it operates. In Section III we describe the fast algorithm, 
and in Section IV we discuss a more sophisticated knapsack cipher. 


ll. AN EXAMPLE OF THE KNAPSACK CIPHER 


A very simple (and insecure) knapsack cipher begins with an “easy” 
knapsack vector generated by a party who wishes to receive encrypted 
data [eq. (29) of Ref. 6], 
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= (1, 2, 4, 8, 17, 35, 68, 142). (1) 


The eight components of E form a super-increasing sequence: Each 
term is larger than the sum of all those preceding it. Let the data to be 
encrypted be represented as a vector with eight binary components, 


= (1, 0, 0, 1, 0, 1, 1, 1). (2) 
To encrypt D using E, form the dot product, 
Sze = E-D = 254. (3) 


The number Sz is an encrypted form of D. 

The super-increasing property of E guarantees that D can be re- 
covered from Sz by subtracting successive components of E (beginning 
with the largest) from Sg and keeping the residue. If a component of 
E is less than or equal to the residue at any stage in the subtraction, 
the corresponding component of D is 1. If a component of E is larger 
than the residue, the corresponding D component is 0 and we try the 
next (smaller) component of E. This process is illustrated below for 
Sz = 254 [eq. (3)], 


a 

1<l<l 9<9 112 oi 

ol ee” = \s nat’ ca 5 \ =68 g \ =e 

0 xX X 112 

D= (1 O O 1 0 1 1 1). 


Of course, E cannot be used for secure encryption, because if E were 
obtained by an eavesdropper he could use it to decrypt any transmitted 
message. The knapsack cipher provides security by transforming E 
into a “hard” knapsack vector H (the public key), which can be used 
for encryption, but which is useless for decryption. To generate this 
transformation, the receiver chooses two secret integers M and W such 
that: (i) M is larger than the sum of all the components in E, and (iz) 
W and M are relatively prime. (This condition means that W is 
invertible modulo M: W~'. W = 1 mod M.) Following Ref. 6, we choose 
M = 291 and W = 176 (which implies W~' = 167). H is generated from 
E by 


H; = W-E; mod M, (4) 
yielding 
= (176; 61, 122, 244, 82, 49, 37, 257). 


In the ideal case H looks like a random sequence; the super-increasing 
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structure of the original E is completely obliterated. H, the public key, 
is sent to the transmitter and need not be kept secret. 
To encrypt D using H, form the dot product as before, 


Su = H-D = 763. (5) 


Su is the encrypted data. If the number of components in H is large, 
say 100 or more, then an eavesdropper, even though he has H and Sz, 
cannot recover D in a reasonable time. The legitimate receiver, how- 
ever, can recover D easily by using the inverse transformation, 


Sz = Sy-W' mod M. (6) 


That is, by using his secret M and W™’, the receiver can convert Sy 
into the number Sz [eq. (3)], the same number that would have been 
obtained if D had been encrypted with E instead of H. Once he has Sz, 
the receiver simply subtracts off successive components of E to recover 
D. 


lll. A FAST DECRYPTION ALGORITHM 


The most time-consuming step of the knapsack cipher is the modular 
multiplication of eq. (6). In practice, the quantities Sy, W~', and M 
might be 100 to 200 bits long, making computation of Sg very slow. 
The calculation can be expedited by considering the n-bit binary 
expansion of Sz, 


Si = bn-102"* + oes + by 2°. (7) 
Substituting eq. (7) into eq. (6), we have 
Sz = [bn-1(2" 'W7 mod M) 
+ +++ + bo(2°W' mod M)]mod M. (8) 


Each term in the square brackets is the product of a binary digit (0 or 
1) and a fixed quantity (in parentheses), which can be computed ahead 
of time and stored in a memory. Evaluation of Sz thus reduces to a 
sequence of table lookups and accumulations, one lookup for each bit 
in Sy. After all the bits in Sy have been processed, the final reduction 
mod M is accomplished by an easy long division. [The division is 
“easy” because each term in eq. (8) can be no bigger than M, so the 
final sum can be no bigger than nM; division by M can therefore be 
accomplished with only approximately logen substract-and-shift oper- 
ations in binary arithmetic. | 

Table I shows the contents of the lookup table required for decryp- 
tion of the example in Section II, along with the binary representation 
of Sy = 763. The value of the sum within the square brackets of eq. 
(8) is seen to be 1127, which is equivalent to 254 in mod 291 arithmetic, 
as required. 
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Table !—Lookup table 
for decryption of the 
example in Section Il 

2*.167 


mod 291 bi 
9 241 1 
8 266 0 
7 133 1 
6 212 1 
5 106 1 
4 53 1 
3 172 1 
2 86 0 
1 43 1 
0 167 1 


Figure 1 shows a block diagram of the decryption process. The basic 
steps of lookup, accumulation, reduction mod M and successive sub- 
traction are pipelined, and within each step most of the processing can 
be performed on all bits in parallel. This architecture results in very 
fast operation, the speed limitation being either the memory access 
time or the accumulator add time, whichever is greater. With Schottky 
TTL and carry-lookahead addition, these times are both in the neigh- 
borhood of 50 ns, so a throughput rate of 10 Mbit/s is reasonable. 

Implementation of the decryption algorithm using very large scale 
integration appears attractive. Most of the circuitry is simply a large 
lookup table, as shown at the top of Fig. 1. Its capacity is determined 
primarily by the number of components in E and the allowed range 
(number of possible values) for each component of E. We can achieve 
reasonable security by using 100 and 2’, respectively, for these two 
parameters; this leads to a value for the modulus M in the neighbor- 
hood of 27”. Since each component of H is less than M, Su [eq. (5)] 
will be less than 27°’. The lookup table must therefore contain 207 
words, each 200 bits long, implying a memory size of approximately 41 
kilobits. Additional memory (~15 kilobits) is required to store the 
components of E. Thus, approximately 56 kilobits of memory and 
some simple arithmetic logic to perform the steps of accumulation, 
long division, and successive subtraction are adequate to implement 
the decryption process. This level of complexity is within the range of 
current vLsI technology.’® 

Finally, we remark that a straightforward implementation of Fig. 1 
may not be the best approach; several modifications of the basic 
decryption algorithm must be investigated. For example, the lookup 
table can be eliminated by calculating the numbers 2*. W~’ mod M 
one-by-one as they are needed for each incoming bit of Sy. Starting 
with W, successive numbers can be generated by a simple left shift 
(and subtraction of M if necessary). 
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IV. ITERATED KNAPSACK TRANSFORMATIONS 


The security of the knapsack cipher is enhanced if iterated (multiple) 
knapsack transformations are employed.’ For example, the “hard” 
vector H [eq. (4)] can be kept secret and used to generate a “harder” 
public vector H’, 


H; = W’-H;mod M’. (9) 
Data can be encrypted with H’ in the usual fashion, 
Sy = H’-D. (10) 


If M’ is chosen to be greater than the sum of all the components of H, 
then data encrypted using H’ may be decrypted using two successive 
inverse transformations having the form of eq. (6). The cost of this 
double-iteration technique in terms of the bandwidth efficiency of the 
cipher is modest. For a 100-component knapsack, the modulus M’ will 





ENCRYPTED 
DATA INPUT 
(n-BIT BLOCK) 





COMPONENTS SUCCESSIVE 
OFE SUBTRACTIONS 


DECRYPTED DATA OUTPUT 





Fig. 1—The fast knapsack decryption algorithm. (Wide arrows signify parallel data 
transfer.) Pipeline architecture and parallel processing contribute to a high throughput 
rate. Hardware implementation would require approximately 56 kilobits of memory and 
a small amount of arithmetic logic. 
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be roughly 100 times bigger than M; thus Sy: will require only about 
seven more bits than Sy would have required. 

We illustrate the double-iteration technique by continuing with the 
example of Section II. Let M’ = 2001 and W’ = 1984, giving (W’)7' = 
1177. From eq. (9) we have 


H’ = (1010, 964, 1928, 1855, 607, 1168, 1372, 1634). (11) 


Encrypting D [eq. (2)] with H’ yields Sy: = 7039. 
Decryption requires two inverse transformations, 


Se = Sy-W 7’ mod M 
= [Sy--(W’)"* mod M’]. W~' mod M 
= [7039-1177 mod 2001]-167 mod 291 
= 763-167 mod 291 
= 254. (12) 


The cascaded inverse transformations in eq. (12) can be executed in 
tandem using the algorithm of Section III. Thus, the decryption 
process will entail a longer total delay (compared with the single- 
iteration case), but the net throughput rate will be essentially un- 
changed. 

We mentioned earlier that straightforward use of the multiple iter- 
ation technique reduces the bandwidth efficiency of the cipher. It is 
possible, however, that for a given level of security, multiple interations 
may actually be more efficient than a single knapsack transformation. 
This is because the enhanced security associated with repeated trans- 
formations might permit a smaller range for the components of E, and 
hence smaller values for the moduli M, M’, etc. The consequent 
reduction in the encrypted block length could offset the seven-bit 
increase normally associated with each iteration. 


V. CONCLUSIONS 


The existence of a fast algorithm for decryption of the knapsack 
cipher means that the advantages of public-key cryptosystems can be 
realized even in high-speed applications. Full integration of the de- 
cryption process onto a single chip appears feasible with current VLSI 
technology. The relationships among cipher security, bandwidth effi- 
ciency, and number of iterations need further investigation. 
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